Introduce BlockSparseMatrixData

A number of algorithms like the SchurEliminator do not need
access to the full BlockSparseMatrix interface. They only
need read only access to the values array and the block structure.

This change introduces, BlockSparseDataMatrix a struct that carries
these two bits of information and modifies the Schur type algorithms
to use it.

What this change will allow us to do, in a subsequent CL is to
take the values array of a BlockSparseMatrix and pair it with
a different blocks structure for subset preconditioning.

Change-Id: I1808f12531b586c9ff4d6a70b3d390c7b0d9f441
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index 4bdbee8..d0c255d 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -133,6 +133,31 @@
   std::unique_ptr<CompressedRowBlockStructure> block_structure_;
 };
 
+// A number of algorithms like the SchurEliminator do not need
+// access to the full BlockSparseMatrix interface. They only
+// need read only access to the values array and the block structure.
+//
+// BlockSparseDataMatrix a struct that carries these two bits of
+// information
+class BlockSparseMatrixData {
+ public:
+  BlockSparseMatrixData(const BlockSparseMatrix& m)
+      : block_structure_(m.block_structure()), values_(m.values()){};
+
+  BlockSparseMatrixData(const CompressedRowBlockStructure* block_structure,
+                        const double* values)
+      : block_structure_(block_structure), values_(values) {}
+
+  const CompressedRowBlockStructure* block_structure() const {
+    return block_structure_;
+  }
+  const double* values() const { return values_; }
+
+ private:
+  const CompressedRowBlockStructure* block_structure_;
+  const double* values_;
+};
+
 }  // namespace internal
 }  // namespace ceres
 
diff --git a/internal/ceres/implicit_schur_complement_test.cc b/internal/ceres/implicit_schur_complement_test.cc
index fb07a52..4e3a598 100644
--- a/internal/ceres/implicit_schur_complement_test.cc
+++ b/internal/ceres/implicit_schur_complement_test.cc
@@ -98,7 +98,8 @@
     lhs->resize(num_schur_rows, num_schur_rows);
     rhs->resize(num_schur_rows);
 
-    eliminator->Eliminate(A_.get(), b_.get(), D, &blhs, rhs->data());
+    eliminator->Eliminate(
+        BlockSparseMatrixData(*A_), b_.get(), D, &blhs, rhs->data());
 
     MatrixRef lhs_ref(blhs.mutable_values(), num_schur_rows, num_schur_rows);
 
@@ -114,7 +115,7 @@
     VectorRef schur_solution(solution->data() + num_cols_ - num_schur_rows,
                              num_schur_rows);
     schur_solution = lhs->selfadjointView<Eigen::Upper>().llt().solve(*rhs);
-    eliminator->BackSubstitute(A_.get(), b_.get(), D,
+    eliminator->BackSubstitute(BlockSparseMatrixData(*A_), b_.get(), D,
                                schur_solution.data(), solution->data());
   }
 
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 27a44a7..0083300 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -156,7 +156,11 @@
   std::fill(x, x + A->num_cols(), 0.0);
   event_logger.AddEvent("Setup");
 
-  eliminator_->Eliminate(A, b, per_solve_options.D, lhs_.get(), rhs_.get());
+  eliminator_->Eliminate(BlockSparseMatrixData(*A),
+                         b,
+                         per_solve_options.D,
+                         lhs_.get(),
+                         rhs_.get());
   event_logger.AddEvent("Eliminate");
 
   double* reduced_solution = x + A->num_cols() - lhs_->num_cols();
@@ -165,7 +169,8 @@
   event_logger.AddEvent("ReducedSolve");
 
   if (summary.termination_type == LINEAR_SOLVER_SUCCESS) {
-    eliminator_->BackSubstitute(A, b, per_solve_options.D, reduced_solution, x);
+    eliminator_->BackSubstitute(
+        BlockSparseMatrixData(*A), b, per_solve_options.D, reduced_solution, x);
     event_logger.AddEvent("BackSubstitute");
   }
 
diff --git a/internal/ceres/schur_eliminator.h b/internal/ceres/schur_eliminator.h
index 7dfc160..66fcb4d 100644
--- a/internal/ceres/schur_eliminator.h
+++ b/internal/ceres/schur_eliminator.h
@@ -195,7 +195,7 @@
   //
   // Since the Schur complement is a symmetric matrix, only the upper
   // triangular part of the Schur complement is computed.
-  virtual void Eliminate(const BlockSparseMatrix* A,
+  virtual void Eliminate(const BlockSparseMatrixData& A,
                          const double* b,
                          const double* D,
                          BlockRandomAccessMatrix* lhs,
@@ -204,7 +204,7 @@
   // Given values for the variables z in the F block of A, solve for
   // the optimal values of the variables y corresponding to the E
   // block in A.
-  virtual void BackSubstitute(const BlockSparseMatrix* A,
+  virtual void BackSubstitute(const BlockSparseMatrixData& A,
                               const double* b,
                               const double* D,
                               const double* z,
@@ -235,12 +235,12 @@
   void Init(int num_eliminate_blocks,
             bool assume_full_rank_ete,
             const CompressedRowBlockStructure* bs) final;
-  void Eliminate(const BlockSparseMatrix* A,
+  void Eliminate(const BlockSparseMatrixData& A,
                  const double* b,
                  const double* D,
                  BlockRandomAccessMatrix* lhs,
                  double* rhs) final;
-  void BackSubstitute(const BlockSparseMatrix* A,
+  void BackSubstitute(const BlockSparseMatrixData& A,
                       const double* b,
                       const double* D,
                       const double* z,
@@ -282,7 +282,7 @@
 
   void ChunkDiagonalBlockAndGradient(
       const Chunk& chunk,
-      const BlockSparseMatrix* A,
+      const BlockSparseMatrixData& A,
       const double* b,
       int row_block_counter,
       typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix* eet,
@@ -291,7 +291,7 @@
       BlockRandomAccessMatrix* lhs);
 
   void UpdateRhs(const Chunk& chunk,
-                 const BlockSparseMatrix* A,
+                 const BlockSparseMatrixData& A,
                  const double* b,
                  int row_block_counter,
                  const double* inverse_ete_g,
@@ -303,17 +303,17 @@
                          const double* buffer,
                          const BufferLayoutType& buffer_layout,
                          BlockRandomAccessMatrix* lhs);
-  void EBlockRowOuterProduct(const BlockSparseMatrix* A,
+  void EBlockRowOuterProduct(const BlockSparseMatrixData& A,
                              int row_block_index,
                              BlockRandomAccessMatrix* lhs);
 
-  void NoEBlockRowsUpdate(const BlockSparseMatrix* A,
+  void NoEBlockRowsUpdate(const BlockSparseMatrixData& A,
                           const double* b,
                           int row_block_counter,
                           BlockRandomAccessMatrix* lhs,
                           double* rhs);
 
-  void NoEBlockRowOuterProduct(const BlockSparseMatrix* A,
+  void NoEBlockRowOuterProduct(const BlockSparseMatrixData& A,
                                int row_block_index,
                                BlockRandomAccessMatrix* lhs);
 
@@ -425,7 +425,7 @@
         e_t_e_inverse_matrices_.begin(), e_t_e_inverse_matrices_.end(), 0.0);
   }
 
-  void Eliminate(const BlockSparseMatrix* A,
+  void Eliminate(const BlockSparseMatrixData& A,
                  const double* b,
                  const double* D,
                  BlockRandomAccessMatrix* lhs_bram,
@@ -442,8 +442,8 @@
     lhs.setZero();
     rhs.setZero();
 
-    const CompressedRowBlockStructure* bs = A->block_structure();
-    const double* values = A->values();
+    const CompressedRowBlockStructure* bs = A.block_structure();
+    const double* values = A.values();
 
     // Add the diagonal to the schur complement.
     if (D != nullptr) {
@@ -566,14 +566,14 @@
   // before this. SchurComplementSolver always does this.
   //
   // y_i = e_t_e_inverse * sum_i e_i^T * (b_i - f_i * z);
-  void BackSubstitute(const BlockSparseMatrix* A,
+  void BackSubstitute(const BlockSparseMatrixData& A,
                       const double* b,
                       const double* D,
                       const double* z_ptr,
                       double* y) override {
     typename EigenTypes<kFBlockSize>::ConstVectorRef z(z_ptr, kFBlockSize);
-    const CompressedRowBlockStructure* bs = A->block_structure();
-    const double* values = A->values();
+    const CompressedRowBlockStructure* bs = A.block_structure();
+    const double* values = A.values();
     Eigen::Matrix<double, kEBlockSize, 1> tmp;
     for (int i = 0; i < chunks_.size(); ++i) {
       const Chunk& chunk = chunks_[i];
diff --git a/internal/ceres/schur_eliminator_benchmark.cc b/internal/ceres/schur_eliminator_benchmark.cc
index 9dcc77f..6307025 100644
--- a/internal/ceres/schur_eliminator_benchmark.cc
+++ b/internal/ceres/schur_eliminator_benchmark.cc
@@ -142,7 +142,7 @@
 
   eliminator->Init(num_e_blocks, true, data.matrix().block_structure());
   for (auto _ : state) {
-    eliminator->Eliminate(&data.matrix(),
+    eliminator->Eliminate(BlockSparseMatrixData(data.matrix()),
                           data.b().data(),
                           data.diagonal().data(),
                           data.mutable_lhs(),
@@ -164,13 +164,13 @@
       SchurEliminatorBase::Create(linear_solver_options));
 
   eliminator->Init(num_e_blocks, true, data.matrix().block_structure());
-  eliminator->Eliminate(&data.matrix(),
+  eliminator->Eliminate(BlockSparseMatrixData(data.matrix()),
                         data.b().data(),
                         data.diagonal().data(),
                         data.mutable_lhs(),
                         data.mutable_rhs()->data());
   for (auto _ : state) {
-    eliminator->BackSubstitute(&data.matrix(),
+    eliminator->BackSubstitute(BlockSparseMatrixData(data.matrix()),
                                data.b().data(),
                                data.diagonal().data(),
                                data.mutable_z()->data(),
@@ -184,7 +184,7 @@
   SchurEliminatorForOneFBlock<2, 3, 6> eliminator;
   eliminator.Init(num_e_blocks, true, data.matrix().block_structure());
   for (auto _ : state) {
-    eliminator.Eliminate(&data.matrix(),
+    eliminator.Eliminate(BlockSparseMatrixData(data.matrix()),
                          data.b().data(),
                          data.diagonal().data(),
                          data.mutable_lhs(),
@@ -197,13 +197,13 @@
   BenchmarkData data(num_e_blocks);
   SchurEliminatorForOneFBlock<2, 3, 6> eliminator;
   eliminator.Init(num_e_blocks, true, data.matrix().block_structure());
-  eliminator.Eliminate(&data.matrix(),
+  eliminator.Eliminate(BlockSparseMatrixData(data.matrix()),
                        data.b().data(),
                        data.diagonal().data(),
                        data.mutable_lhs(),
                        data.mutable_rhs()->data());
   for (auto _ : state) {
-    eliminator.BackSubstitute(&data.matrix(),
+    eliminator.BackSubstitute(BlockSparseMatrixData(data.matrix()),
                               data.b().data(),
                               data.diagonal().data(),
                               data.mutable_z()->data(),
diff --git a/internal/ceres/schur_eliminator_impl.h b/internal/ceres/schur_eliminator_impl.h
index 5d65869..bd0881e 100644
--- a/internal/ceres/schur_eliminator_impl.h
+++ b/internal/ceres/schur_eliminator_impl.h
@@ -104,6 +104,13 @@
     lhs_num_rows += bs->cols[i].size;
   }
 
+  // TODO(sameeragarwal): Now that we may have subset block structure,
+  // we need to make sure that we account for the fact that somep
+  // point blocks only have a "diagonal" row and nothing more.
+  //
+  // This likely requires a slightly different algorithm, which works
+  // off of the number of elimination blocks.
+
   int r = 0;
   // Iterate over the row blocks of A, and detect the chunks. The
   // matrix should already have been ordered so that all rows
@@ -145,7 +152,7 @@
       ++chunk.size;
     }
 
-    CHECK_GT(chunk.size, 0);
+    CHECK_GT(chunk.size, 0); // This check will need to be resolved.
     r += chunk.size;
   }
   const Chunk& chunk = chunks_.back();
@@ -169,7 +176,7 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
-Eliminate(const BlockSparseMatrix* A,
+Eliminate(const BlockSparseMatrixData& A,
           const double* b,
           const double* D,
           BlockRandomAccessMatrix* lhs,
@@ -181,7 +188,7 @@
     }
   }
 
-  const CompressedRowBlockStructure* bs = A->block_structure();
+  const CompressedRowBlockStructure* bs = A.block_structure();
   const int num_col_blocks = bs->cols.size();
 
   // Add the diagonal to the schur complement.
@@ -303,12 +310,13 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
-BackSubstitute(const BlockSparseMatrix* A,
+BackSubstitute(const BlockSparseMatrixData& A,
                const double* b,
                const double* D,
                const double* z,
                double* y) {
-  const CompressedRowBlockStructure* bs = A->block_structure();
+  const CompressedRowBlockStructure* bs = A.block_structure();
+  const double* values = A.values();
 
   ParallelFor(
       context_,
@@ -333,7 +341,6 @@
       ete.setZero();
     }
 
-    const double* values = A->values();
     for (int j = 0; j < chunk.size; ++j) {
       const CompressedRow& row = bs->rows[chunk.start + j];
       const Cell& e_cell = row.cells.front();
@@ -380,17 +387,17 @@
 void
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
 UpdateRhs(const Chunk& chunk,
-          const BlockSparseMatrix* A,
+          const BlockSparseMatrixData& A,
           const double* b,
           int row_block_counter,
           const double* inverse_ete_g,
           double* rhs) {
-  const CompressedRowBlockStructure* bs = A->block_structure();
+  const CompressedRowBlockStructure* bs = A.block_structure();
+  const double* values = A.values();
+
   const int e_block_id = bs->rows[chunk.start].cells.front().block_id;
   const int e_block_size = bs->cols[e_block_id].size;
-
   int b_pos = bs->rows[row_block_counter].block.position;
-  const double* values = A->values();
   for (int j = 0; j < chunk.size; ++j) {
     const CompressedRow& row = bs->rows[row_block_counter + j];
     const Cell& e_cell = row.cells.front();
@@ -441,14 +448,15 @@
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
 ChunkDiagonalBlockAndGradient(
     const Chunk& chunk,
-    const BlockSparseMatrix* A,
+    const BlockSparseMatrixData& A,
     const double* b,
     int row_block_counter,
     typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix* ete,
     double* g,
     double* buffer,
     BlockRandomAccessMatrix* lhs) {
-  const CompressedRowBlockStructure* bs = A->block_structure();
+  const CompressedRowBlockStructure* bs = A.block_structure();
+  const double* values = A.values();
 
   int b_pos = bs->rows[row_block_counter].block.position;
   const int e_block_size = ete->rows();
@@ -457,7 +465,6 @@
   // contribution of its F blocks to the Schur complement, the
   // contribution of its E block to the matrix EE' (ete), and the
   // corresponding block in the gradient vector.
-  const double* values = A->values();
   for (int j = 0; j < chunk.size; ++j) {
     const CompressedRow& row = bs->rows[row_block_counter + j];
 
@@ -558,13 +565,13 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
-NoEBlockRowsUpdate(const BlockSparseMatrix* A,
+NoEBlockRowsUpdate(const BlockSparseMatrixData& A,
                    const double* b,
                    int row_block_counter,
                    BlockRandomAccessMatrix* lhs,
                    double* rhs) {
-  const CompressedRowBlockStructure* bs = A->block_structure();
-  const double* values = A->values();
+  const CompressedRowBlockStructure* bs = A.block_structure();
+  const double* values = A.values();
   for (; row_block_counter < bs->rows.size(); ++row_block_counter) {
     NoEBlockRowOuterProduct(A, row_block_counter, lhs);
     if (!rhs) {
@@ -601,12 +608,13 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
-NoEBlockRowOuterProduct(const BlockSparseMatrix* A,
+NoEBlockRowOuterProduct(const BlockSparseMatrixData& A,
                         int row_block_index,
                         BlockRandomAccessMatrix* lhs) {
-  const CompressedRowBlockStructure* bs = A->block_structure();
+  const CompressedRowBlockStructure* bs = A.block_structure();
+  const double* values = A.values();
+
   const CompressedRow& row = bs->rows[row_block_index];
-  const double* values = A->values();
   for (int i = 0; i < row.cells.size(); ++i) {
     const int block1 = row.cells[i].block_id - num_eliminate_blocks_;
     DCHECK_GE(block1, 0);
@@ -654,12 +662,13 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void
 SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
-EBlockRowOuterProduct(const BlockSparseMatrix* A,
+EBlockRowOuterProduct(const BlockSparseMatrixData& A,
                       int row_block_index,
                       BlockRandomAccessMatrix* lhs) {
-  const CompressedRowBlockStructure* bs = A->block_structure();
+  const CompressedRowBlockStructure* bs = A.block_structure();
+  const double* values = A.values();
+
   const CompressedRow& row = bs->rows[row_block_index];
-  const double* values = A->values();
   for (int i = 1; i < row.cells.size(); ++i) {
     const int block1 = row.cells[i].block_id - num_eliminate_blocks_;
     DCHECK_GE(block1, 0);
diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc
index 645b25a..6383ced 100644
--- a/internal/ceres/schur_eliminator_test.cc
+++ b/internal/ceres/schur_eliminator_test.cc
@@ -57,8 +57,8 @@
 class SchurEliminatorTest : public ::testing::Test {
  protected:
   void SetUpFromId(int id) {
-    std::unique_ptr<LinearLeastSquaresProblem>
-        problem(CreateLinearLeastSquaresProblemFromId(id));
+    std::unique_ptr<LinearLeastSquaresProblem> problem(
+        CreateLinearLeastSquaresProblemFromId(id));
     CHECK(problem != nullptr);
     SetupHelper(problem.get());
   }
@@ -85,7 +85,7 @@
     A->ToDenseMatrix(&J);
     VectorRef f(b.get(), J.rows());
 
-    Matrix H  =  (D.cwiseProduct(D)).asDiagonal();
+    Matrix H = (D.cwiseProduct(D)).asDiagonal();
     H.noalias() += J.transpose() * J;
 
     const Vector g = J.transpose() * f;
@@ -101,28 +101,21 @@
     sol_expected.setZero();
 
     Matrix P = H.block(0, 0, num_eliminate_cols, num_eliminate_cols);
-    Matrix Q = H.block(0,
-                       num_eliminate_cols,
-                       num_eliminate_cols,
-                       schur_size);
-    Matrix R = H.block(num_eliminate_cols,
-                       num_eliminate_cols,
-                       schur_size,
-                       schur_size);
+    Matrix Q = H.block(0, num_eliminate_cols, num_eliminate_cols, schur_size);
+    Matrix R =
+        H.block(num_eliminate_cols, num_eliminate_cols, schur_size, schur_size);
     int row = 0;
     const CompressedRowBlockStructure* bs = A->block_structure();
     for (int i = 0; i < num_eliminate_blocks; ++i) {
-      const int block_size =  bs->cols[i].size;
-      P.block(row, row,  block_size, block_size) =
-          P
-          .block(row, row,  block_size, block_size)
-          .llt()
-          .solve(Matrix::Identity(block_size, block_size));
+      const int block_size = bs->cols[i].size;
+      P.block(row, row, block_size, block_size) =
+          P.block(row, row, block_size, block_size)
+              .llt()
+              .solve(Matrix::Identity(block_size, block_size));
       row += block_size;
     }
 
-    lhs_expected
-        .triangularView<Eigen::Upper>() = R - Q.transpose() * P * Q;
+    lhs_expected.triangularView<Eigen::Upper>() = R - Q.transpose() * P * Q;
     rhs_expected =
         g.tail(schur_size) - Q.transpose() * P * g.head(num_eliminate_cols);
     sol_expected = H.llt().solve(g);
@@ -161,20 +154,18 @@
     eliminator.reset(SchurEliminatorBase::Create(options));
     const bool kFullRankETE = true;
     eliminator->Init(num_eliminate_blocks, kFullRankETE, A->block_structure());
-    eliminator->Eliminate(A.get(), b.get(), diagonal.data(), &lhs, rhs.data());
+    eliminator->Eliminate(
+        BlockSparseMatrixData(*A), b.get(), diagonal.data(), &lhs, rhs.data());
 
     MatrixRef lhs_ref(lhs.mutable_values(), lhs.num_rows(), lhs.num_cols());
-    Vector reduced_sol  =
-        lhs_ref
-        .selfadjointView<Eigen::Upper>()
-        .llt()
-        .solve(rhs);
+    Vector reduced_sol =
+        lhs_ref.selfadjointView<Eigen::Upper>().llt().solve(rhs);
 
     // Solution to the linear least squares problem.
     Vector sol(num_cols);
     sol.setZero();
     sol.tail(schur_size) = reduced_sol;
-    eliminator->BackSubstitute(A.get(),
+    eliminator->BackSubstitute(BlockSparseMatrixData(*A),
                                b.get(),
                                diagonal.data(),
                                reduced_sol.data(),
@@ -183,9 +174,11 @@
     Matrix delta = (lhs_ref - lhs_expected).selfadjointView<Eigen::Upper>();
     double diff = delta.norm();
     EXPECT_NEAR(diff / lhs_expected.norm(), 0.0, relative_tolerance);
-    EXPECT_NEAR((rhs - rhs_expected).norm() / rhs_expected.norm(), 0.0,
+    EXPECT_NEAR((rhs - rhs_expected).norm() / rhs_expected.norm(),
+                0.0,
                 relative_tolerance);
-    EXPECT_NEAR((sol - sol_expected).norm() / sol_expected.norm(), 0.0,
+    EXPECT_NEAR((sol - sol_expected).norm() / sol_expected.norm(),
+                0.0,
                 relative_tolerance);
   }
 
@@ -324,39 +317,53 @@
     std::unique_ptr<SchurEliminatorBase> eliminator(
         SchurEliminatorBase::Create(linear_solver_options));
     eliminator->Init(num_e_blocks, true, matrix.block_structure());
-    eliminator->Eliminate(&matrix, b.data(), diagonal.data(), &expected_lhs,
+    eliminator->Eliminate(BlockSparseMatrixData(matrix),
+                          b.data(),
+                          diagonal.data(),
+                          &expected_lhs,
                           expected_rhs.data());
-    eliminator->BackSubstitute(&matrix, b.data(), diagonal.data(), f_sol.data(),
+    eliminator->BackSubstitute(BlockSparseMatrixData(matrix),
+                               b.data(),
+                               diagonal.data(),
+                               f_sol.data(),
                                actual_e_sol.data());
   }
 
   {
     SchurEliminatorForOneFBlock<2, 3, 6> eliminator;
     eliminator.Init(num_e_blocks, true, matrix.block_structure());
-    eliminator.Eliminate(&matrix, b.data(), diagonal.data(), &actual_lhs,
+    eliminator.Eliminate(BlockSparseMatrixData(matrix),
+                         b.data(),
+                         diagonal.data(),
+                         &actual_lhs,
                          actual_rhs.data());
-    eliminator.BackSubstitute(&matrix, b.data(), diagonal.data(), f_sol.data(),
+    eliminator.BackSubstitute(BlockSparseMatrixData(matrix),
+                              b.data(),
+                              diagonal.data(),
+                              f_sol.data(),
                               expected_e_sol.data());
   }
-  ConstMatrixRef actual_lhsref(actual_lhs.values(), actual_lhs.num_cols(),
-                               actual_lhs.num_cols());
-  ConstMatrixRef expected_lhsref(expected_lhs.values(), actual_lhs.num_cols(),
-                                 actual_lhs.num_cols());
+  ConstMatrixRef actual_lhsref(
+      actual_lhs.values(), actual_lhs.num_cols(), actual_lhs.num_cols());
+  ConstMatrixRef expected_lhsref(
+      expected_lhs.values(), actual_lhs.num_cols(), actual_lhs.num_cols());
 
   EXPECT_NEAR((actual_lhsref - expected_lhsref).norm() / expected_lhsref.norm(),
-              0.0, 1e-12)
+              0.0,
+              1e-12)
       << "expected: \n"
       << expected_lhsref << "\nactual: \n"
       << actual_lhsref;
 
-  EXPECT_NEAR((actual_rhs - expected_rhs).norm() / expected_rhs.norm(), 0.0,
-              1e-12)
+  EXPECT_NEAR(
+      (actual_rhs - expected_rhs).norm() / expected_rhs.norm(), 0.0, 1e-12)
       << "expected: \n"
       << expected_rhs << "\nactual: \n"
       << actual_rhs;
 
   EXPECT_NEAR((actual_e_sol - expected_e_sol).norm() / expected_e_sol.norm(),
-              0.0, 1e-12)
+              0.0,
+              1e-12)
       << "expected: \n"
       << expected_e_sol << "\nactual: \n"
       << actual_e_sol;
diff --git a/internal/ceres/schur_jacobi_preconditioner.cc b/internal/ceres/schur_jacobi_preconditioner.cc
index 1500650..89d770b 100644
--- a/internal/ceres/schur_jacobi_preconditioner.cc
+++ b/internal/ceres/schur_jacobi_preconditioner.cc
@@ -49,9 +49,8 @@
   CHECK_GT(options_.elimination_groups.size(), 1);
   CHECK_GT(options_.elimination_groups[0], 0);
   const int num_blocks = bs.cols.size() - options_.elimination_groups[0];
-  CHECK_GT(num_blocks, 0)
-      << "Jacobian should have at least 1 f_block for "
-      << "SCHUR_JACOBI preconditioner.";
+  CHECK_GT(num_blocks, 0) << "Jacobian should have at least 1 f_block for "
+                          << "SCHUR_JACOBI preconditioner.";
   CHECK(options_.context != NULL);
 
   std::vector<int> blocks(num_blocks);
@@ -63,8 +62,7 @@
   InitEliminator(bs);
 }
 
-SchurJacobiPreconditioner::~SchurJacobiPreconditioner() {
-}
+SchurJacobiPreconditioner::~SchurJacobiPreconditioner() {}
 
 // Initialize the SchurEliminator.
 void SchurJacobiPreconditioner::InitEliminator(
@@ -89,7 +87,8 @@
   CHECK_GT(num_rows, 0);
 
   // Compute a subset of the entries of the Schur complement.
-  eliminator_->Eliminate(&A, nullptr, D, m_.get(), nullptr);
+  eliminator_->Eliminate(
+      BlockSparseMatrixData(A), nullptr, D, m_.get(), nullptr);
   m_->Invert();
   return true;
 }
@@ -99,9 +98,7 @@
   m_->RightMultiply(x, y);
 }
 
-int SchurJacobiPreconditioner::num_rows() const {
-  return m_->num_rows();
-}
+int SchurJacobiPreconditioner::num_rows() const { return m_->num_rows(); }
 
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index ed4afb6..2b60fac 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -185,7 +185,7 @@
 // The cluster_membership_ vector is updated to indicate cluster
 // memberships for each camera block.
 void VisibilityBasedPreconditioner::ClusterCameras(
-    const vector<set<int> >& visibility) {
+    const vector<set<int>>& visibility) {
   std::unique_ptr<WeightedGraph<int>> schur_complement_graph(
       CreateSchurComplementGraph(visibility));
   CHECK(schur_complement_graph != nullptr);
@@ -342,7 +342,8 @@
   CHECK_GT(num_rows, 0);
 
   // Compute a subset of the entries of the Schur complement.
-  eliminator_->Eliminate(&A, nullptr, D, m_.get(), nullptr);
+  eliminator_->Eliminate(
+      BlockSparseMatrixData(A), nullptr, D, m_.get(), nullptr);
 
   // Try factorizing the matrix. For CLUSTER_JACOBI, this should
   // always succeed modulo some numerical/conditioning problems. For
@@ -464,7 +465,7 @@
 // each vertex.
 void VisibilityBasedPreconditioner::ForestToClusterPairs(
     const WeightedGraph<int>& forest,
-    std::unordered_set<pair<int, int>, pair_hash >* cluster_pairs) const {
+    std::unordered_set<pair<int, int>, pair_hash>* cluster_pairs) const {
   CHECK(cluster_pairs != nullptr);
   cluster_pairs->clear();
   const std::unordered_set<int>& vertices = forest.vertices();