diff --git a/internal/ceres/evaluation_benchmark.cc b/internal/ceres/evaluation_benchmark.cc
index 9849b5a..9717d0f 100644
--- a/internal/ceres/evaluation_benchmark.cc
+++ b/internal/ceres/evaluation_benchmark.cc
@@ -147,19 +147,28 @@
 
   const PartitionedView* PartitionedMatrixViewJacobian(
       const LinearSolver::Options& options) {
-    auto block_sparse = BlockSparseJacobian(options.context);
+    auto block_sparse = BlockSparseJacobianWithTranspose(options.context);
     partitioned_view_jacobian =
         std::make_unique<PartitionedView>(options, *block_sparse);
     return partitioned_view_jacobian.get();
   }
 
-  const PartitionedView* PartitionedMatrixViewJacobianWithTranspose(
-      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();
+  BlockSparseMatrix* BlockDiagonalEtE(const LinearSolver::Options& options) {
+    if (!block_diagonal_ete) {
+      auto partitioned_view =
+          PartitionedMatrixViewJacobian(options);
+      block_diagonal_ete = partitioned_view->CreateBlockDiagonalEtE();
+    }
+    return block_diagonal_ete.get();
+  }
+
+  BlockSparseMatrix* BlockDiagonalFtF(const LinearSolver::Options& options) {
+    if (!block_diagonal_ftf) {
+      auto partitioned_view =
+          PartitionedMatrixViewJacobian(options);
+      block_diagonal_ftf = partitioned_view->CreateBlockDiagonalFtF();
+    }
+    return block_diagonal_ftf.get();
   }
 
   Vector parameters;
@@ -169,7 +178,8 @@
   std::unique_ptr<BlockSparseMatrix> block_sparse_jacobian_with_transpose;
   std::unique_ptr<CompressedRowSparseMatrix> crs_jacobian;
   std::unique_ptr<PartitionedView> partitioned_view_jacobian;
-  std::unique_ptr<PartitionedView> partitioned_view_jacobian_with_transpose;
+  std::unique_ptr<BlockSparseMatrix> block_diagonal_ete;
+  std::unique_ptr<BlockSparseMatrix> block_diagonal_ftf;
 };
 
 static void Residuals(benchmark::State& state,
@@ -262,7 +272,7 @@
   options.num_threads = state.range(0);
   options.elimination_groups.push_back(data->bal_problem->num_points());
   options.context = context;
-  auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(options);
+  auto jacobian = data->PartitionedMatrixViewJacobian(options);
 
   Vector y = Vector::Zero(jacobian->num_cols_f());
   Vector x = Vector::Random(jacobian->num_rows());
@@ -298,7 +308,7 @@
   options.num_threads = state.range(0);
   options.elimination_groups.push_back(data->bal_problem->num_points());
   options.context = context;
-  auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(options);
+  auto jacobian = data->PartitionedMatrixViewJacobian(options);
 
   Vector y = Vector::Zero(jacobian->num_cols_e());
   Vector x = Vector::Random(jacobian->num_rows());
@@ -309,6 +319,36 @@
   CHECK_GT(y.squaredNorm(), 0.);
 }
 
+static void PMVUpdateBlockDiagonalEtE(benchmark::State& state,
+                                      BALData* data,
+                                      ContextImpl* 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);
+  auto block_diagonal_ete = data->BlockDiagonalEtE(options);
+
+  for (auto _ : state) {
+    jacobian->UpdateBlockDiagonalEtE(block_diagonal_ete);
+  }
+}
+
+static void PMVUpdateBlockDiagonalFtF(benchmark::State& state,
+                                      BALData* data,
+                                      ContextImpl* 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);
+  auto block_diagonal_ftf = data->BlockDiagonalFtF(options);
+
+  for (auto _ : state) {
+    jacobian->UpdateBlockDiagonalFtF(block_diagonal_ftf);
+  }
+}
+
 static void JacobianRightMultiplyAndAccumulate(benchmark::State& state,
                                                BALData* data,
                                                ContextImpl* context) {
@@ -485,6 +525,30 @@
         ->Arg(8)
         ->Arg(16);
 
+    const std::string name_update_block_diagonal_ftf =
+        "PMVUpdateBlockDiagonalFtF<" + path + ">";
+    ::benchmark::RegisterBenchmark(name_update_block_diagonal_ftf.c_str(),
+                                   ceres::internal::PMVUpdateBlockDiagonalFtF,
+                                   data,
+                                   &context)
+        ->Arg(1)
+        ->Arg(2)
+        ->Arg(4)
+        ->Arg(8)
+        ->Arg(16);
+
+    const std::string name_update_block_diagonal_ete =
+        "PMVUpdateBlockDiagonalEtE<" + path + ">";
+    ::benchmark::RegisterBenchmark(name_update_block_diagonal_ete.c_str(),
+                                   ceres::internal::PMVUpdateBlockDiagonalEtE,
+                                   data,
+                                   &context)
+        ->Arg(1)
+        ->Arg(2)
+        ->Arg(4)
+        ->Arg(8)
+        ->Arg(16);
+
 #ifndef CERES_NO_CUDA
     const std::string name_right_product_cuda =
         "JacobianRightMultiplyAndAccumulateCuda<" + path + ">";
diff --git a/internal/ceres/partitioned_matrix_view.h b/internal/ceres/partitioned_matrix_view.h
index f34fd4a..380b576 100644
--- a/internal/ceres/partitioned_matrix_view.h
+++ b/internal/ceres/partitioned_matrix_view.h
@@ -166,7 +166,15 @@
   std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalEtE() const final;
   std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalFtF() const final;
   void UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const final;
+  void UpdateBlockDiagonalEtESingleThreaded(
+      BlockSparseMatrix* block_diagonal) const;
+  void UpdateBlockDiagonalEtEMultiThreaded(
+      BlockSparseMatrix* block_diagonal) const;
   void UpdateBlockDiagonalFtF(BlockSparseMatrix* block_diagonal) const final;
+  void UpdateBlockDiagonalFtFSingleThreaded(
+      BlockSparseMatrix* block_diagonal) const;
+  void UpdateBlockDiagonalFtFMultiThreaded(
+      BlockSparseMatrix* block_diagonal) const;
   // clang-format off
   int num_col_blocks_e() const final { return num_col_blocks_e_;  }
   int num_col_blocks_f() const final { return num_col_blocks_f_;  }
diff --git a/internal/ceres/partitioned_matrix_view_impl.h b/internal/ceres/partitioned_matrix_view_impl.h
index 1ca2e14..e685f5b 100644
--- a/internal/ceres/partitioned_matrix_view_impl.h
+++ b/internal/ceres/partitioned_matrix_view_impl.h
@@ -444,10 +444,10 @@
 //
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const {
-  const CompressedRowBlockStructure* bs = matrix_.block_structure();
-  const CompressedRowBlockStructure* block_diagonal_structure =
-      block_diagonal->block_structure();
+    UpdateBlockDiagonalEtESingleThreaded(
+        BlockSparseMatrix* block_diagonal) const {
+  auto bs = matrix_.block_structure();
+  auto block_diagonal_structure = block_diagonal->block_structure();
 
   block_diagonal->SetZero();
   const double* values = matrix_.values();
@@ -470,6 +470,57 @@
   }
 }
 
+template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
+void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
+    UpdateBlockDiagonalEtEMultiThreaded(
+        BlockSparseMatrix* block_diagonal) const {
+  auto transpose_block_structure = matrix_.transpose_block_structure();
+  CHECK(transpose_block_structure != nullptr);
+  auto block_diagonal_structure = block_diagonal->block_structure();
+
+  const double* values = matrix_.values();
+  double* values_diagonal = block_diagonal->mutable_values();
+  ParallelFor(
+      options_.context,
+      0,
+      num_col_blocks_e_,
+      options_.num_threads,
+      [values,
+       transpose_block_structure,
+       values_diagonal,
+       block_diagonal_structure](int col_block_id) {
+        int cell_position =
+            block_diagonal_structure->rows[col_block_id].cells[0].position;
+        double* cell_values = values_diagonal + cell_position;
+        int col_block_size =
+            transpose_block_structure->rows[col_block_id].block.size;
+        auto& cells = transpose_block_structure->rows[col_block_id].cells;
+        MatrixRef(cell_values, col_block_size, col_block_size).setZero();
+
+        for (auto& c : cells) {
+          int row_block_size = transpose_block_structure->cols[c.block_id].size;
+          // clang-format off
+          MatrixTransposeMatrixMultiply<kRowBlockSize, kEBlockSize, kRowBlockSize, kEBlockSize, 1>(
+            values + c.position, row_block_size, col_block_size,
+            values + c.position, row_block_size, col_block_size,
+            cell_values, 0, 0, col_block_size, col_block_size);
+          // clang-format on
+        }
+      },
+      e_cols_partition_);
+}
+
+template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
+void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
+    UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const {
+  if (options_.num_threads == 1) {
+    UpdateBlockDiagonalEtESingleThreaded(block_diagonal);
+  } else {
+    CHECK(options_.context != nullptr);
+    UpdateBlockDiagonalEtEMultiThreaded(block_diagonal);
+  }
+}
+
 // Similar to the code in RightMultiplyAndAccumulateF, except instead of the
 // matrix vector multiply its an outer product.
 //
@@ -477,10 +528,10 @@
 //
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    UpdateBlockDiagonalFtF(BlockSparseMatrix* block_diagonal) const {
-  const CompressedRowBlockStructure* bs = matrix_.block_structure();
-  const CompressedRowBlockStructure* block_diagonal_structure =
-      block_diagonal->block_structure();
+    UpdateBlockDiagonalFtFSingleThreaded(
+        BlockSparseMatrix* block_diagonal) const {
+  auto bs = matrix_.block_structure();
+  auto block_diagonal_structure = block_diagonal->block_structure();
 
   block_diagonal->SetZero();
   const double* values = matrix_.values();
@@ -527,4 +578,82 @@
   }
 }
 
+template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
+void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
+    UpdateBlockDiagonalFtFMultiThreaded(
+        BlockSparseMatrix* block_diagonal) const {
+  auto transpose_block_structure = matrix_.transpose_block_structure();
+  CHECK(transpose_block_structure != nullptr);
+  auto block_diagonal_structure = block_diagonal->block_structure();
+
+  const double* values = matrix_.values();
+  double* values_diagonal = block_diagonal->mutable_values();
+
+  const int num_col_blocks_e = num_col_blocks_e_;
+  const int num_row_blocks_e = num_row_blocks_e_;
+  ParallelFor(
+      options_.context,
+      num_col_blocks_e_,
+      num_col_blocks_e + num_col_blocks_f_,
+      options_.num_threads,
+      [transpose_block_structure,
+       block_diagonal_structure,
+       num_col_blocks_e,
+       num_row_blocks_e,
+       values,
+       values_diagonal](int col_block_id) {
+        const int col_block_size =
+            transpose_block_structure->rows[col_block_id].block.size;
+        const int diagonal_block_id = col_block_id - num_col_blocks_e;
+        const int cell_position =
+            block_diagonal_structure->rows[diagonal_block_id].cells[0].position;
+        double* cell_values = values_diagonal + cell_position;
+
+        MatrixRef(cell_values, col_block_size, col_block_size).setZero();
+
+        auto& cells = transpose_block_structure->rows[col_block_id].cells;
+        const int num_cells = cells.size();
+        int i = 0;
+        for (; i < num_cells; ++i) {
+          auto& cell = cells[i];
+          const int row_block_id = cell.block_id;
+          if (row_block_id >= num_row_blocks_e) break;
+          const int row_block_size =
+              transpose_block_structure->cols[row_block_id].size;
+          // clang-format off
+          MatrixTransposeMatrixMultiply
+              <kRowBlockSize, kFBlockSize, kRowBlockSize, kFBlockSize, 1>(
+                  values + cell.position, row_block_size, col_block_size,
+                  values + cell.position, row_block_size, col_block_size,
+                  cell_values, 0, 0, col_block_size, col_block_size);
+          // clang-format on
+        }
+        for (; i < num_cells; ++i) {
+          auto& cell = cells[i];
+          const int row_block_id = cell.block_id;
+          const int row_block_size =
+              transpose_block_structure->cols[row_block_id].size;
+          // clang-format off
+          MatrixTransposeMatrixMultiply
+              <Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic, 1>(
+                  values + cell.position, row_block_size, col_block_size,
+                  values + cell.position, row_block_size, col_block_size,
+                  cell_values, 0, 0, col_block_size, col_block_size);
+          // clang-format on
+        }
+      },
+      f_cols_partition_);
+}
+
+template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
+void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
+    UpdateBlockDiagonalFtF(BlockSparseMatrix* block_diagonal) const {
+  if (options_.num_threads == 1) {
+    UpdateBlockDiagonalFtFSingleThreaded(block_diagonal);
+  } else {
+    CHECK(options_.context != nullptr);
+    UpdateBlockDiagonalFtFMultiThreaded(block_diagonal);
+  }
+}
+
 }  // namespace ceres::internal
diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc
index 3d10b1c..f74351e 100644
--- a/internal/ceres/partitioned_matrix_view_test.cc
+++ b/internal/ceres/partitioned_matrix_view_test.cc
@@ -49,63 +49,6 @@
 
 const double kEpsilon = 1e-14;
 
-class PartitionedMatrixViewTest : public ::testing::Test {
- protected:
-  void SetUp() final {
-    std::unique_ptr<LinearLeastSquaresProblem> problem =
-        CreateLinearLeastSquaresProblemFromId(2);
-    CHECK(problem != nullptr);
-    A_ = std::move(problem->A);
-
-    num_cols_ = A_->num_cols();
-    num_rows_ = A_->num_rows();
-    num_eliminate_blocks_ = problem->num_eliminate_blocks;
-    LinearSolver::Options options;
-    options.elimination_groups.push_back(num_eliminate_blocks_);
-    pmv_ = PartitionedMatrixViewBase::Create(
-        options, *down_cast<BlockSparseMatrix*>(A_.get()));
-  }
-
-  double RandDouble() { return distribution_(prng_); }
-
-  int num_rows_;
-  int num_cols_;
-  int num_eliminate_blocks_;
-  std::unique_ptr<SparseMatrix> A_;
-  std::unique_ptr<PartitionedMatrixViewBase> pmv_;
-  std::mt19937 prng_;
-  std::uniform_real_distribution<double> distribution_ =
-      std::uniform_real_distribution<double>(0.0, 1.0);
-};
-
-TEST_F(PartitionedMatrixViewTest, BlockDiagonalEtE) {
-  std::unique_ptr<BlockSparseMatrix> block_diagonal_ee(
-      pmv_->CreateBlockDiagonalEtE());
-  const CompressedRowBlockStructure* bs = block_diagonal_ee->block_structure();
-
-  EXPECT_EQ(block_diagonal_ee->num_rows(), 2);
-  EXPECT_EQ(block_diagonal_ee->num_cols(), 2);
-  EXPECT_EQ(bs->cols.size(), 2);
-  EXPECT_EQ(bs->rows.size(), 2);
-
-  EXPECT_NEAR(block_diagonal_ee->values()[0], 10.0, kEpsilon);
-  EXPECT_NEAR(block_diagonal_ee->values()[1], 155.0, kEpsilon);
-}
-
-TEST_F(PartitionedMatrixViewTest, BlockDiagonalFtF) {
-  std::unique_ptr<BlockSparseMatrix> block_diagonal_ff(
-      pmv_->CreateBlockDiagonalFtF());
-  const CompressedRowBlockStructure* bs = block_diagonal_ff->block_structure();
-
-  EXPECT_EQ(block_diagonal_ff->num_rows(), 3);
-  EXPECT_EQ(block_diagonal_ff->num_cols(), 3);
-  EXPECT_EQ(bs->cols.size(), 3);
-  EXPECT_EQ(bs->rows.size(), 3);
-  EXPECT_NEAR(block_diagonal_ff->values()[0], 70.0, kEpsilon);
-  EXPECT_NEAR(block_diagonal_ff->values()[1], 17.0, kEpsilon);
-  EXPECT_NEAR(block_diagonal_ff->values()[2], 37.0, kEpsilon);
-}
-
 // Param = <problem_id, num_threads>
 using Param = ::testing::tuple<int, int>;
 
@@ -116,7 +59,7 @@
   return ss.str();
 }
 
-class PartitionedMatrixViewSpMVTest : public ::testing::TestWithParam<Param> {
+class PartitionedMatrixViewTest : public ::testing::TestWithParam<Param> {
  protected:
   void SetUp() final {
     const int problem_id = ::testing::get<0>(GetParam());
@@ -132,6 +75,10 @@
     options_.elimination_groups.push_back(problem->num_eliminate_blocks);
     pmv_ = PartitionedMatrixViewBase::Create(options_, *block_sparse);
 
+    LinearSolver::Options options_single_threaded = options_;
+    options_single_threaded.num_threads = 1;
+    pmv_single_threaded_ = 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() -
@@ -147,12 +94,13 @@
   std::unique_ptr<LinearLeastSquaresProblem> problem_;
   std::unique_ptr<SparseMatrix> A_;
   std::unique_ptr<PartitionedMatrixViewBase> pmv_;
+  std::unique_ptr<PartitionedMatrixViewBase> pmv_single_threaded_;
   std::mt19937 prng_;
   std::uniform_real_distribution<double> distribution_ =
       std::uniform_real_distribution<double>(0.0, 1.0);
 };
 
-TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateE) {
+TEST_P(PartitionedMatrixViewTest, RightMultiplyAndAccumulateE) {
   Vector x1(pmv_->num_cols_e());
   Vector x2(pmv_->num_cols());
   x2.setZero();
@@ -172,7 +120,7 @@
   }
 }
 
-TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateF) {
+TEST_P(PartitionedMatrixViewTest, RightMultiplyAndAccumulateF) {
   Vector x1(pmv_->num_cols_f());
   Vector x2(pmv_->num_cols());
   x2.setZero();
@@ -192,7 +140,7 @@
   }
 }
 
-TEST_P(PartitionedMatrixViewSpMVTest, LeftMultiplyAndAccumulate) {
+TEST_P(PartitionedMatrixViewTest, LeftMultiplyAndAccumulate) {
   Vector x = Vector::Zero(pmv_->num_rows());
   for (int i = 0; i < pmv_->num_rows(); ++i) {
     x(i) = RandDouble();
@@ -215,9 +163,109 @@
   }
 }
 
+TEST_P(PartitionedMatrixViewTest, BlockDiagonalFtF) {
+  std::unique_ptr<BlockSparseMatrix> block_diagonal_ff(
+      pmv_->CreateBlockDiagonalFtF());
+  const auto bs_diagonal = block_diagonal_ff->block_structure();
+  const int num_rows = pmv_->num_rows();
+  const int num_cols_f = pmv_->num_cols_f();
+  const int num_cols_e = pmv_->num_cols_e();
+  const int num_col_blocks_f = pmv_->num_col_blocks_f();
+  const int num_col_blocks_e = pmv_->num_col_blocks_e();
+
+  CHECK_EQ(block_diagonal_ff->num_rows(), num_cols_f);
+  CHECK_EQ(block_diagonal_ff->num_cols(), num_cols_f);
+
+  EXPECT_EQ(bs_diagonal->cols.size(), num_col_blocks_f);
+  EXPECT_EQ(bs_diagonal->rows.size(), num_col_blocks_f);
+
+  Matrix EF;
+  A_->ToDenseMatrix(&EF);
+  const auto F = EF.topRightCorner(num_rows, num_cols_f);
+
+  Matrix expected_FtF = F.transpose() * F;
+  Matrix actual_FtF;
+  block_diagonal_ff->ToDenseMatrix(&actual_FtF);
+
+  // FtF might be not block-diagonal
+  auto bs = down_cast<BlockSparseMatrix*>(A_.get())->block_structure();
+  for (int i = 0; i < num_col_blocks_f; ++i) {
+    const auto col_block_f = bs->cols[num_col_blocks_e + i];
+    const int block_size = col_block_f.size;
+    const int block_pos = col_block_f.position - num_cols_e;
+    const auto cell_expected =
+        expected_FtF.block(block_pos, block_pos, block_size, block_size);
+    auto cell_actual =
+        actual_FtF.block(block_pos, block_pos, block_size, block_size);
+    cell_actual -= cell_expected;
+    EXPECT_NEAR(cell_actual.norm(), 0., kEpsilon);
+  }
+  // There should be nothing remaining outside block-diagonal
+  EXPECT_NEAR(actual_FtF.norm(), 0., kEpsilon);
+}
+
+TEST_P(PartitionedMatrixViewTest, BlockDiagonalEtE) {
+  std::unique_ptr<BlockSparseMatrix> block_diagonal_ee(
+      pmv_->CreateBlockDiagonalEtE());
+  const CompressedRowBlockStructure* bs = block_diagonal_ee->block_structure();
+  const int num_rows = pmv_->num_rows();
+  const int num_cols_e = pmv_->num_cols_e();
+  const int num_col_blocks_e = pmv_->num_col_blocks_e();
+
+  CHECK_EQ(block_diagonal_ee->num_rows(), num_cols_e);
+  CHECK_EQ(block_diagonal_ee->num_cols(), num_cols_e);
+
+  EXPECT_EQ(bs->cols.size(), num_col_blocks_e);
+  EXPECT_EQ(bs->rows.size(), num_col_blocks_e);
+
+  Matrix EF;
+  A_->ToDenseMatrix(&EF);
+  const auto E = EF.topLeftCorner(num_rows, num_cols_e);
+
+  Matrix expected_EtE = E.transpose() * E;
+  Matrix actual_EtE;
+  block_diagonal_ee->ToDenseMatrix(&actual_EtE);
+
+  EXPECT_NEAR((expected_EtE - actual_EtE).norm(), 0., kEpsilon);
+}
+
+TEST_P(PartitionedMatrixViewTest, UpdateBlockDiagonalEtE) {
+  std::unique_ptr<BlockSparseMatrix> block_diagonal_ete(
+      pmv_->CreateBlockDiagonalEtE());
+  const CompressedRowBlockStructure* bs = block_diagonal_ete->block_structure();
+  const int num_cols = pmv_->num_cols_e();
+
+  Matrix multi_threaded(num_cols, num_cols);
+  pmv_->UpdateBlockDiagonalEtE(block_diagonal_ete.get());
+  block_diagonal_ete->ToDenseMatrix(&multi_threaded);
+
+  Matrix single_threaded(num_cols, num_cols);
+  pmv_single_threaded_->UpdateBlockDiagonalEtE(block_diagonal_ete.get());
+  block_diagonal_ete->ToDenseMatrix(&single_threaded);
+
+  EXPECT_NEAR((multi_threaded - single_threaded).norm(), 0., kEpsilon);
+}
+
+TEST_P(PartitionedMatrixViewTest, UpdateBlockDiagonalFtF) {
+  std::unique_ptr<BlockSparseMatrix> block_diagonal_ftf(
+      pmv_->CreateBlockDiagonalFtF());
+  const CompressedRowBlockStructure* bs = block_diagonal_ftf->block_structure();
+  const int num_cols = pmv_->num_cols_f();
+
+  Matrix multi_threaded(num_cols, num_cols);
+  pmv_->UpdateBlockDiagonalFtF(block_diagonal_ftf.get());
+  block_diagonal_ftf->ToDenseMatrix(&multi_threaded);
+
+  Matrix single_threaded(num_cols, num_cols);
+  pmv_single_threaded_->UpdateBlockDiagonalFtF(block_diagonal_ftf.get());
+  block_diagonal_ftf->ToDenseMatrix(&single_threaded);
+
+  EXPECT_NEAR((multi_threaded - single_threaded).norm(), 0., kEpsilon);
+}
+
 INSTANTIATE_TEST_SUITE_P(
     ParallelProducts,
-    PartitionedMatrixViewSpMVTest,
+    PartitionedMatrixViewTest,
     ::testing::Combine(::testing::Values(2, 4, 6),
                        ::testing::Values(1, 2, 3, 4, 5, 6, 7, 8)),
     ParamInfoToString);