Unify Block handling across matrix types

Previously some matrices used Block to keep track of
row/column block sizes and some would just use ints, and
then compute the position of each row and column from it.

By uniformly using Block everywhere, we reduce duplicate
computation and data copies.

I also cleaned up a bunch of c++17 related stuff as I edited
these files.

Change-Id: I4c86b1593fd4c91f9057fbb38314f62f303e0477
diff --git a/internal/ceres/block_jacobi_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc
index 7aaa605..db70abf 100644
--- a/internal/ceres/block_jacobi_preconditioner.cc
+++ b/internal/ceres/block_jacobi_preconditioner.cc
@@ -42,13 +42,8 @@
 
 BlockSparseJacobiPreconditioner::BlockSparseJacobiPreconditioner(
     const BlockSparseMatrix& A) {
-  const CompressedRowBlockStructure* bs = A.block_structure();
-  std::vector<int> blocks(bs->cols.size());
-  for (int i = 0; i < blocks.size(); ++i) {
-    blocks[i] = bs->cols[i].size;
-  }
-
-  m_ = std::make_unique<BlockRandomAccessDiagonalMatrix>(blocks);
+  m_ = std::make_unique<BlockRandomAccessDiagonalMatrix>(
+      A.block_structure()->cols);
 }
 
 BlockSparseJacobiPreconditioner::~BlockSparseJacobiPreconditioner() = default;
@@ -104,11 +99,7 @@
 
   // Compute the number of non-zeros in the preconditioner. This is needed so
   // that we can construct the CompressedRowSparseMatrix.
-  int m_nnz = 0;
-  for (int col_block_size : col_blocks) {
-    m_nnz += col_block_size * col_block_size;
-  }
-
+  const int m_nnz = SumSquaredSizes(col_blocks);
   m_ = std::make_unique<CompressedRowSparseMatrix>(
       A.num_cols(), A.num_cols(), m_nnz);
 
@@ -123,14 +114,13 @@
     // Not that the because of the way the CompressedRowSparseMatrix format
     // works, the entire diagonal block is laid out contiguously in memory as a
     // row-major matrix. We will use this when updating the block.
-    const int col_block_size = col_blocks[i];
-    for (int j = 0; j < col_block_size; ++j) {
-      for (int k = 0; k < col_block_size; ++k, ++idx) {
-        m_cols[idx] = col + k;
+    auto& block = col_blocks[i];
+    for (int j = 0; j < block.size; ++j) {
+      for (int k = 0; k < block.size; ++k, ++idx) {
+        m_cols[idx] = block.position + k;
       }
-      m_rows[col + j + 1] = idx;
+      m_rows[block.position + j + 1] = idx;
     }
-    col += col_block_size;
   }
 
   CHECK_EQ(m_rows[A.num_cols()], m_nnz);
@@ -154,15 +144,14 @@
   const int* m_rows = m_->rows();
 
   const int num_rows = A.num_rows();
-  int r = 0;
   for (int i = 0; i < num_row_blocks; ++i) {
-    const int row_block_size = row_blocks[i];
-    const int row_nnz = a_rows[r + 1] - a_rows[r];
-    ConstMatrixRef row_block(a_values + a_rows[r], row_block_size, row_nnz);
-    int idx = a_rows[r];
+    const int row = row_blocks[i].position;
+    const int row_block_size = row_blocks[i].size;
+    const int row_nnz = a_rows[row + 1] - a_rows[row];
+    ConstMatrixRef row_block(a_values + a_rows[row], row_block_size, row_nnz);
     int c = 0;
     while (c < row_nnz) {
-      const int idx = a_rows[r] + c;
+      const int idx = a_rows[row] + c;
       const int col = a_cols[idx];
       const int col_block_size = m_rows[col + 1] - m_rows[col];
 
@@ -177,11 +166,11 @@
       m.noalias() += b.transpose() * b;
       c += col_block_size;
     }
-    r += row_block_size;
   }
 
-  for (int i = 0, col = 0; i < num_col_blocks; ++i) {
-    const int col_block_size = m_rows[col + 1] - m_rows[col];
+  for (int i = 0; i < num_col_blocks; ++i) {
+    const int col = col_blocks[i].position;
+    const int col_block_size = col_blocks[i].size;
     MatrixRef m(m_values + m_rows[col], col_block_size, col_block_size);
 
     if (D != nullptr) {
@@ -189,8 +178,9 @@
           ConstVectorRef(D + col, col_block_size).array().square().matrix();
     }
 
+    // TODO(sameeragarwal): Deal with Cholesky inversion failure here and
+    // elsewhere.
     m = m.llt().solve(Matrix::Identity(col_block_size, col_block_size));
-    col += col_block_size;
   }
 
   return true;
diff --git a/internal/ceres/block_jacobi_preconditioner_test.cc b/internal/ceres/block_jacobi_preconditioner_test.cc
index db60c7b..fa051e4 100644
--- a/internal/ceres/block_jacobi_preconditioner_test.cc
+++ b/internal/ceres/block_jacobi_preconditioner_test.cc
@@ -123,7 +123,7 @@
 
     const auto& col_blocks = jacobian->col_blocks();
     for (int i = 0, col = 0; i < col_blocks.size(); ++i) {
-      const int block_size = col_blocks[i];
+      const int block_size = col_blocks[i].size;
       int idx = m.rows()[col];
       for (int j = 0; j < block_size; ++j) {
         EXPECT_EQ(m.rows()[col + j + 1] - m.rows()[col + j], block_size);
diff --git a/internal/ceres/block_random_access_dense_matrix.cc b/internal/ceres/block_random_access_dense_matrix.cc
index 341d7fe..78c2146 100644
--- a/internal/ceres/block_random_access_dense_matrix.cc
+++ b/internal/ceres/block_random_access_dense_matrix.cc
@@ -38,17 +38,11 @@
 namespace ceres::internal {
 
 BlockRandomAccessDenseMatrix::BlockRandomAccessDenseMatrix(
-    const std::vector<int>& blocks) {
-  const int num_blocks = blocks.size();
-  block_layout_.resize(num_blocks, 0);
-  num_rows_ = 0;
-  for (int i = 0; i < num_blocks; ++i) {
-    block_layout_[i] = num_rows_;
-    num_rows_ += blocks[i];
-  }
-
+    std::vector<Block> blocks)
+    : blocks_(std::move(blocks)) {
+  const int num_blocks = blocks_.size();
+  num_rows_ = NumScalarEntries(blocks_);
   values_ = std::make_unique<double[]>(num_rows_ * num_rows_);
-
   cell_infos_ = std::make_unique<CellInfo[]>(num_blocks * num_blocks);
   for (int i = 0; i < num_blocks * num_blocks; ++i) {
     cell_infos_[i].values = values_.get();
@@ -67,11 +61,11 @@
                                                 int* col,
                                                 int* row_stride,
                                                 int* col_stride) {
-  *row = block_layout_[row_block_id];
-  *col = block_layout_[col_block_id];
+  *row = blocks_[row_block_id].position;
+  *col = blocks_[col_block_id].position;
   *row_stride = num_rows_;
   *col_stride = num_rows_;
-  return &cell_infos_[row_block_id * block_layout_.size() + col_block_id];
+  return &cell_infos_[row_block_id * blocks_.size() + col_block_id];
 }
 
 // Assume that the user does not hold any locks on any cell blocks
diff --git a/internal/ceres/block_random_access_dense_matrix.h b/internal/ceres/block_random_access_dense_matrix.h
index c8c066a..085a54d 100644
--- a/internal/ceres/block_random_access_dense_matrix.h
+++ b/internal/ceres/block_random_access_dense_matrix.h
@@ -35,6 +35,7 @@
 #include <vector>
 
 #include "ceres/block_random_access_matrix.h"
+#include "ceres/block_structure.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
 
@@ -55,7 +56,7 @@
  public:
   // blocks is a vector of block sizes. The resulting matrix has
   // blocks.size() * blocks.size() cells.
-  explicit BlockRandomAccessDenseMatrix(const std::vector<int>& blocks);
+  explicit BlockRandomAccessDenseMatrix(std::vector<Block> blocks);
   BlockRandomAccessDenseMatrix(const BlockRandomAccessDenseMatrix&) = delete;
   void operator=(const BlockRandomAccessDenseMatrix&) = delete;
 
@@ -86,7 +87,7 @@
 
  private:
   int num_rows_;
-  std::vector<int> block_layout_;
+  std::vector<Block> blocks_;
   std::unique_ptr<double[]> values_;
   std::unique_ptr<CellInfo[]> cell_infos_;
 };
diff --git a/internal/ceres/block_random_access_dense_matrix_test.cc b/internal/ceres/block_random_access_dense_matrix_test.cc
index c85c388..b0cd5b9 100644
--- a/internal/ceres/block_random_access_dense_matrix_test.cc
+++ b/internal/ceres/block_random_access_dense_matrix_test.cc
@@ -35,23 +35,22 @@
 #include "ceres/internal/eigen.h"
 #include "gtest/gtest.h"
 
-namespace ceres {
-namespace internal {
+namespace ceres::internal {
 
 TEST(BlockRandomAccessDenseMatrix, GetCell) {
-  std::vector<int> blocks;
-  blocks.push_back(3);
-  blocks.push_back(4);
-  blocks.push_back(5);
+  std::vector<Block> blocks;
+  blocks.emplace_back(3, 0);
+  blocks.emplace_back(4, 3);
+  blocks.emplace_back(5, 7);
   const int num_rows = 3 + 4 + 5;
   BlockRandomAccessDenseMatrix m(blocks);
   EXPECT_EQ(m.num_rows(), num_rows);
   EXPECT_EQ(m.num_cols(), num_rows);
 
-  int row_idx = 0;
   for (int i = 0; i < blocks.size(); ++i) {
-    int col_idx = 0;
+    const int row_idx = blocks[i].position;
     for (int j = 0; j < blocks.size(); ++j) {
+      const int col_idx = blocks[j].position;
       int row;
       int col;
       int row_stride;
@@ -63,17 +62,15 @@
       EXPECT_EQ(col, col_idx);
       EXPECT_EQ(row_stride, 3 + 4 + 5);
       EXPECT_EQ(col_stride, 3 + 4 + 5);
-      col_idx += blocks[j];
     }
-    row_idx += blocks[i];
   }
 }
 
 TEST(BlockRandomAccessDenseMatrix, WriteCell) {
-  std::vector<int> blocks;
-  blocks.push_back(3);
-  blocks.push_back(4);
-  blocks.push_back(5);
+  std::vector<Block> blocks;
+  blocks.emplace_back(3, 0);
+  blocks.emplace_back(4, 3);
+  blocks.emplace_back(5, 7);
   const int num_rows = 3 + 4 + 5;
 
   BlockRandomAccessDenseMatrix m(blocks);
@@ -87,29 +84,26 @@
       int col_stride;
       CellInfo* cell = m.GetCell(i, j, &row, &col, &row_stride, &col_stride);
       MatrixRef(cell->values, row_stride, col_stride)
-          .block(row, col, blocks[i], blocks[j]) =
-          (i + 1) * (j + 1) * Matrix::Ones(blocks[i], blocks[j]);
+          .block(row, col, blocks[i].size, blocks[j].size) =
+          (i + 1) * (j + 1) * Matrix::Ones(blocks[i].size, blocks[j].size);
     }
   }
 
   // Check the values in the array are correct by going over the
   // entries of each block manually.
-  int row_idx = 0;
   for (int i = 0; i < blocks.size(); ++i) {
-    int col_idx = 0;
+    const int row_idx = blocks[i].position;
     for (int j = 0; j < blocks.size(); ++j) {
+      const int col_idx = blocks[j].position;
       // Check the values of this block.
-      for (int r = 0; r < blocks[i]; ++r) {
-        for (int c = 0; c < blocks[j]; ++c) {
+      for (int r = 0; r < blocks[i].size; ++r) {
+        for (int c = 0; c < blocks[j].size; ++c) {
           int pos = row_idx * num_rows + col_idx;
           EXPECT_EQ(m.values()[pos], (i + 1) * (j + 1));
         }
       }
-      col_idx += blocks[j];
     }
-    row_idx += blocks[i];
   }
 }
 
-}  // namespace internal
-}  // namespace ceres
+}  // namespace ceres::internal
diff --git a/internal/ceres/block_random_access_diagonal_matrix.cc b/internal/ceres/block_random_access_diagonal_matrix.cc
index 006713f..2e29604 100644
--- a/internal/ceres/block_random_access_diagonal_matrix.cc
+++ b/internal/ceres/block_random_access_diagonal_matrix.cc
@@ -45,24 +45,13 @@
 
 namespace ceres::internal {
 
-using std::vector;
-
 // TODO(sameeragarwal): Drop the dependence on TripletSparseMatrix.
 
 BlockRandomAccessDiagonalMatrix::BlockRandomAccessDiagonalMatrix(
-    const vector<int>& blocks)
-    : blocks_(blocks) {
-  // Build the row/column layout vector and count the number of scalar
-  // rows/columns.
-  int num_cols = 0;
-  int num_nonzeros = 0;
-  vector<int> block_positions;
-  for (int block_size : blocks_) {
-    block_positions.push_back(num_cols);
-    num_cols += block_size;
-    num_nonzeros += block_size * block_size;
-  }
-
+    std::vector<Block> blocks)
+    : blocks_(std::move(blocks)) {
+  const int num_cols = NumScalarEntries(blocks_);
+  const int num_nonzeros = SumSquaredSizes(blocks_);
   VLOG(1) << "Matrix Size [" << num_cols << "," << num_cols << "] "
           << num_nonzeros;
 
@@ -74,14 +63,12 @@
   double* values = tsm_->mutable_values();
 
   int pos = 0;
-  for (int i = 0; i < blocks_.size(); ++i) {
-    const int block_size = blocks_[i];
+  for (auto& block : blocks_) {
     layout_.push_back(new CellInfo(values + pos));
-    const int block_begin = block_positions[i];
-    for (int r = 0; r < block_size; ++r) {
-      for (int c = 0; c < block_size; ++c, ++pos) {
-        rows[pos] = block_begin + r;
-        cols[pos] = block_begin + c;
+    for (int r = 0; r < block.size; ++r) {
+      for (int c = 0; c < block.size; ++c, ++pos) {
+        rows[pos] = block.position + r;
+        cols[pos] = block.position + c;
       }
     }
   }
@@ -102,7 +89,7 @@
   if (row_block_id != col_block_id) {
     return nullptr;
   }
-  const int stride = blocks_[row_block_id];
+  const int stride = blocks_[row_block_id].size;
 
   // Each cell is stored contiguously as its own little dense matrix.
   *row = 0;
@@ -122,11 +109,11 @@
 
 void BlockRandomAccessDiagonalMatrix::Invert() {
   double* values = tsm_->mutable_values();
-  for (int block_size : blocks_) {
-    MatrixRef block(values, block_size, block_size);
-    block = block.selfadjointView<Eigen::Upper>().llt().solve(
-        Matrix::Identity(block_size, block_size));
-    values += block_size * block_size;
+  for (auto& block : blocks_) {
+    MatrixRef b(values, block.size, block.size);
+    b = b.selfadjointView<Eigen::Upper>().llt().solve(
+        Matrix::Identity(block.size, block.size));
+    values += block.size * block.size;
   }
 }
 
@@ -135,12 +122,12 @@
   CHECK(x != nullptr);
   CHECK(y != nullptr);
   const double* values = tsm_->values();
-  for (int block_size : blocks_) {
-    ConstMatrixRef block(values, block_size, block_size);
-    VectorRef(y, block_size).noalias() += block * ConstVectorRef(x, block_size);
-    x += block_size;
-    y += block_size;
-    values += block_size * block_size;
+  for (auto& block : blocks_) {
+    ConstMatrixRef b(values, block.size, block.size);
+    VectorRef(y, block.size).noalias() += b * ConstVectorRef(x, block.size);
+    x += block.size;
+    y += block.size;
+    values += block.size * block.size;
   }
 }
 
diff --git a/internal/ceres/block_random_access_diagonal_matrix.h b/internal/ceres/block_random_access_diagonal_matrix.h
index 2a726a0..d058396 100644
--- a/internal/ceres/block_random_access_diagonal_matrix.h
+++ b/internal/ceres/block_random_access_diagonal_matrix.h
@@ -37,6 +37,7 @@
 #include <vector>
 
 #include "ceres/block_random_access_matrix.h"
+#include "ceres/block_structure.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
 #include "ceres/triplet_sparse_matrix.h"
@@ -50,7 +51,7 @@
     : public BlockRandomAccessMatrix {
  public:
   // blocks is an array of block sizes.
-  explicit BlockRandomAccessDiagonalMatrix(const std::vector<int>& blocks);
+  explicit BlockRandomAccessDiagonalMatrix(std::vector<Block> blocks);
   BlockRandomAccessDiagonalMatrix(const BlockRandomAccessDiagonalMatrix&) =
       delete;
   void operator=(const BlockRandomAccessDiagonalMatrix&) = delete;
@@ -86,7 +87,7 @@
 
  private:
   // row/column block sizes.
-  const std::vector<int> blocks_;
+  const std::vector<Block> blocks_;
   std::vector<CellInfo*> layout_;
 
   // The underlying matrix object which actually stores the cells.
diff --git a/internal/ceres/block_random_access_diagonal_matrix_test.cc b/internal/ceres/block_random_access_diagonal_matrix_test.cc
index 37e1f88..5bf47ff 100644
--- a/internal/ceres/block_random_access_diagonal_matrix_test.cc
+++ b/internal/ceres/block_random_access_diagonal_matrix_test.cc
@@ -39,16 +39,16 @@
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
-namespace ceres {
-namespace internal {
+namespace ceres::internal {
 
 class BlockRandomAccessDiagonalMatrixTest : public ::testing::Test {
  public:
   void SetUp() override {
-    std::vector<int> blocks;
-    blocks.push_back(3);
-    blocks.push_back(4);
-    blocks.push_back(5);
+    std::vector<Block> blocks;
+    blocks.emplace_back(3, 0);
+    blocks.emplace_back(4, 3);
+    blocks.emplace_back(5, 7);
+
     const int num_rows = 3 + 4 + 5;
     num_nonzeros_ = 3 * 3 + 4 * 4 + 5 * 5;
 
@@ -78,15 +78,20 @@
         EXPECT_TRUE(cell != nullptr);
         EXPECT_EQ(row, 0);
         EXPECT_EQ(col, 0);
-        EXPECT_EQ(row_stride, blocks[row_block_id]);
-        EXPECT_EQ(col_stride, blocks[col_block_id]);
+        EXPECT_EQ(row_stride, blocks[row_block_id].size);
+        EXPECT_EQ(col_stride, blocks[col_block_id].size);
 
         // Write into the block
         MatrixRef(cell->values, row_stride, col_stride)
-            .block(row, col, blocks[row_block_id], blocks[col_block_id]) =
+            .block(row,
+                   col,
+                   blocks[row_block_id].size,
+                   blocks[col_block_id].size) =
             (row_block_id + 1) * (col_block_id + 1) *
-                Matrix::Ones(blocks[row_block_id], blocks[col_block_id]) +
-            Matrix::Identity(blocks[row_block_id], blocks[row_block_id]);
+                Matrix::Ones(blocks[row_block_id].size,
+                             blocks[col_block_id].size) +
+            Matrix::Identity(blocks[row_block_id].size,
+                             blocks[row_block_id].size);
       }
     }
   }
@@ -160,5 +165,4 @@
   EXPECT_NEAR((expected_inverse - dense).norm(), 0.0, kTolerance);
 }
 
-}  // namespace internal
-}  // namespace ceres
+}  // namespace ceres::internal
diff --git a/internal/ceres/block_random_access_sparse_matrix.cc b/internal/ceres/block_random_access_sparse_matrix.cc
index 2df4c71..f27fc08 100644
--- a/internal/ceres/block_random_access_sparse_matrix.cc
+++ b/internal/ceres/block_random_access_sparse_matrix.cc
@@ -43,32 +43,21 @@
 
 namespace ceres::internal {
 
-using std::make_pair;
-using std::pair;
-using std::set;
-using std::vector;
-
 BlockRandomAccessSparseMatrix::BlockRandomAccessSparseMatrix(
-    const vector<int>& blocks, const set<pair<int, int>>& block_pairs)
+    const std::vector<Block>& blocks,
+    const std::set<std::pair<int, int>>& block_pairs)
     : kMaxRowBlocks(10 * 1000 * 1000), blocks_(blocks) {
   CHECK_LT(blocks.size(), kMaxRowBlocks);
 
-  // Build the row/column layout vector and count the number of scalar
-  // rows/columns.
-  int num_cols = 0;
-  block_positions_.reserve(blocks_.size());
-  for (int block_size : blocks_) {
-    block_positions_.push_back(num_cols);
-    num_cols += block_size;
-  }
+  const int num_cols = NumScalarEntries(blocks);
 
   // Count the number of scalar non-zero entries and build the layout
   // object for looking into the values array of the
   // TripletSparseMatrix.
   int num_nonzeros = 0;
   for (const auto& block_pair : block_pairs) {
-    const int row_block_size = blocks_[block_pair.first];
-    const int col_block_size = blocks_[block_pair.second];
+    const int row_block_size = blocks_[block_pair.first].size;
+    const int col_block_size = blocks_[block_pair.second].size;
     num_nonzeros += row_block_size * col_block_size;
   }
 
@@ -84,8 +73,8 @@
 
   int pos = 0;
   for (const auto& block_pair : block_pairs) {
-    const int row_block_size = blocks_[block_pair.first];
-    const int col_block_size = blocks_[block_pair.second];
+    const int row_block_size = blocks_[block_pair.first].size;
+    const int col_block_size = blocks_[block_pair.second].size;
     cell_values_.emplace_back(block_pair, values + pos);
     layout_[IntPairToLong(block_pair.first, block_pair.second)] =
         new CellInfo(values + pos);
@@ -96,14 +85,14 @@
   for (const auto& block_pair : block_pairs) {
     const int row_block_id = block_pair.first;
     const int col_block_id = block_pair.second;
-    const int row_block_size = blocks_[row_block_id];
-    const int col_block_size = blocks_[col_block_id];
+    const int row_block_size = blocks_[row_block_id].size;
+    const int col_block_size = blocks_[col_block_id].size;
     int pos =
         layout_[IntPairToLong(row_block_id, col_block_id)]->values - values;
     for (int r = 0; r < row_block_size; ++r) {
       for (int c = 0; c < col_block_size; ++c, ++pos) {
-        rows[pos] = block_positions_[row_block_id] + r;
-        cols[pos] = block_positions_[col_block_id] + c;
+        rows[pos] = blocks_[row_block_id].position + r;
+        cols[pos] = blocks_[col_block_id].position + c;
         values[pos] = 1.0;
         DCHECK_LT(rows[pos], tsm_->num_rows());
         DCHECK_LT(cols[pos], tsm_->num_rows());
@@ -134,8 +123,8 @@
   // Each cell is stored contiguously as its own little dense matrix.
   *row = 0;
   *col = 0;
-  *row_stride = blocks_[row_block_id];
-  *col_stride = blocks_[col_block_id];
+  *row_stride = blocks_[row_block_id].size;
+  *col_stride = blocks_[col_block_id].size;
   return it->second;
 }
 
@@ -151,12 +140,12 @@
     const double* x, double* y) const {
   for (const auto& cell_position_and_data : cell_values_) {
     const int row = cell_position_and_data.first.first;
-    const int row_block_size = blocks_[row];
-    const int row_block_pos = block_positions_[row];
+    const int row_block_size = blocks_[row].size;
+    const int row_block_pos = blocks_[row].position;
 
     const int col = cell_position_and_data.first.second;
-    const int col_block_size = blocks_[col];
-    const int col_block_pos = block_positions_[col];
+    const int col_block_size = blocks_[col].size;
+    const int col_block_pos = blocks_[col].position;
 
     MatrixVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
         cell_position_and_data.second,
diff --git a/internal/ceres/block_random_access_sparse_matrix.h b/internal/ceres/block_random_access_sparse_matrix.h
index fe2b13c..7161e99 100644
--- a/internal/ceres/block_random_access_sparse_matrix.h
+++ b/internal/ceres/block_random_access_sparse_matrix.h
@@ -39,6 +39,7 @@
 #include <vector>
 
 #include "ceres/block_random_access_matrix.h"
+#include "ceres/block_structure.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
 #include "ceres/small_blas.h"
@@ -58,7 +59,7 @@
   // <row_block_id, col_block_id> pairs to identify the non-zero cells
   // of this matrix.
   BlockRandomAccessSparseMatrix(
-      const std::vector<int>& blocks,
+      const std::vector<Block>& blocks,
       const std::set<std::pair<int, int>>& block_pairs);
   BlockRandomAccessSparseMatrix(const BlockRandomAccessSparseMatrix&) = delete;
   void operator=(const BlockRandomAccessSparseMatrix&) = delete;
@@ -106,8 +107,7 @@
   const int64_t kMaxRowBlocks;
 
   // row/column block sizes.
-  const std::vector<int> blocks_;
-  std::vector<int> block_positions_;
+  const std::vector<Block> blocks_;
 
   // A mapping from <row_block_id, col_block_id> to the position in
   // the values array of tsm_ where the block is stored.
diff --git a/internal/ceres/block_random_access_sparse_matrix_test.cc b/internal/ceres/block_random_access_sparse_matrix_test.cc
index 605bfbe..78a6c3e 100644
--- a/internal/ceres/block_random_access_sparse_matrix_test.cc
+++ b/internal/ceres/block_random_access_sparse_matrix_test.cc
@@ -32,40 +32,36 @@
 
 #include <limits>
 #include <memory>
+#include <set>
+#include <utility>
 #include <vector>
 
 #include "ceres/internal/eigen.h"
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
-namespace ceres {
-namespace internal {
-
-using std::make_pair;
-using std::pair;
-using std::set;
-using std::vector;
+namespace ceres::internal {
 
 TEST(BlockRandomAccessSparseMatrix, GetCell) {
-  vector<int> blocks;
-  blocks.push_back(3);
-  blocks.push_back(4);
-  blocks.push_back(5);
+  std::vector<Block> blocks;
+  blocks.emplace_back(3, 0);
+  blocks.emplace_back(4, 3);
+  blocks.emplace_back(5, 7);
   const int num_rows = 3 + 4 + 5;
 
-  set<pair<int, int>> block_pairs;
+  std::set<std::pair<int, int>> block_pairs;
   int num_nonzeros = 0;
-  block_pairs.insert(make_pair(0, 0));
-  num_nonzeros += blocks[0] * blocks[0];
+  block_pairs.emplace(0, 0);
+  num_nonzeros += blocks[0].size * blocks[0].size;
 
-  block_pairs.insert(make_pair(1, 1));
-  num_nonzeros += blocks[1] * blocks[1];
+  block_pairs.emplace(1, 1);
+  num_nonzeros += blocks[1].size * blocks[1].size;
 
-  block_pairs.insert(make_pair(1, 2));
-  num_nonzeros += blocks[1] * blocks[2];
+  block_pairs.emplace(1, 2);
+  num_nonzeros += blocks[1].size * blocks[2].size;
 
-  block_pairs.insert(make_pair(0, 2));
-  num_nonzeros += blocks[2] * blocks[0];
+  block_pairs.emplace(0, 2);
+  num_nonzeros += blocks[2].size * blocks[0].size;
 
   BlockRandomAccessSparseMatrix m(blocks, block_pairs);
   EXPECT_EQ(m.num_rows(), num_rows);
@@ -83,14 +79,14 @@
     EXPECT_TRUE(cell != nullptr);
     EXPECT_EQ(row, 0);
     EXPECT_EQ(col, 0);
-    EXPECT_EQ(row_stride, blocks[row_block_id]);
-    EXPECT_EQ(col_stride, blocks[col_block_id]);
+    EXPECT_EQ(row_stride, blocks[row_block_id].size);
+    EXPECT_EQ(col_stride, blocks[col_block_id].size);
 
     // Write into the block
     MatrixRef(cell->values, row_stride, col_stride)
-        .block(row, col, blocks[row_block_id], blocks[col_block_id]) =
+        .block(row, col, blocks[row_block_id].size, blocks[col_block_id].size) =
         (row_block_id + 1) * (col_block_id + 1) *
-        Matrix::Ones(blocks[row_block_id], blocks[col_block_id]);
+        Matrix::Ones(blocks[row_block_id].size, blocks[col_block_id].size);
   }
 
   const TripletSparseMatrix* tsm = m.matrix();
@@ -138,10 +134,10 @@
 class BlockRandomAccessSparseMatrixTest : public ::testing::Test {
  public:
   void SetUp() final {
-    vector<int> blocks;
-    blocks.push_back(1);
-    set<pair<int, int>> block_pairs;
-    block_pairs.insert(make_pair(0, 0));
+    std::vector<Block> blocks;
+    blocks.emplace_back(1, 0);
+    std::set<std::pair<int, int>> block_pairs;
+    block_pairs.emplace(0, 0);
     m_ = std::make_unique<BlockRandomAccessSparseMatrix>(blocks, block_pairs);
   }
 
@@ -179,5 +175,4 @@
   CheckLongToIntPair();
 }
 
-}  // namespace internal
-}  // namespace ceres
+}  // namespace ceres::internal
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index b1447f8..0401696 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -45,8 +45,6 @@
 
 namespace ceres::internal {
 
-using std::vector;
-
 BlockSparseMatrix::BlockSparseMatrix(
     CompressedRowBlockStructure* block_structure)
     : num_rows_(0),
@@ -66,7 +64,7 @@
     int row_block_size = block_structure_->rows[i].block.size;
     num_rows_ += row_block_size;
 
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    const std::vector<Cell>& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -96,7 +94,7 @@
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_pos = block_structure_->rows[i].block.position;
     int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    const std::vector<Cell>& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -119,7 +117,7 @@
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_pos = block_structure_->rows[i].block.position;
     int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    auto& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -139,7 +137,7 @@
   VectorRef(x, num_cols_).setZero();
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    auto& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -156,7 +154,7 @@
 
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    auto& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -181,14 +179,14 @@
   auto& row_blocks = *crs_matrix->mutable_row_blocks();
   row_blocks.resize(num_row_blocks);
   for (int i = 0; i < num_row_blocks; ++i) {
-    row_blocks[i] = block_structure_->rows[i].block.size;
+    row_blocks[i] = block_structure_->rows[i].block;
   }
 
   int num_col_blocks = block_structure_->cols.size();
   auto& col_blocks = *crs_matrix->mutable_col_blocks();
   col_blocks.resize(num_col_blocks);
   for (int i = 0; i < num_col_blocks; ++i) {
-    col_blocks[i] = block_structure_->cols[i].size;
+    col_blocks[i] = block_structure_->cols[i];
   }
 }
 
@@ -202,7 +200,7 @@
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_pos = block_structure_->rows[i].block.position;
     int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    auto& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -225,7 +223,7 @@
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     int row_block_pos = block_structure_->rows[i].block.position;
     int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    const auto& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       int col_block_id = cell.block_id;
       int col_block_size = block_structure_->cols[col_block_id].size;
@@ -254,7 +252,7 @@
   for (int i = 0; i < block_structure_->rows.size(); ++i) {
     const int row_block_pos = block_structure_->rows[i].block.position;
     const int row_block_size = block_structure_->rows[i].block.size;
-    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    const auto& cells = block_structure_->rows[i].cells;
     for (const auto& cell : cells) {
       const int col_block_id = cell.block_id;
       const int col_block_size = block_structure_->cols[col_block_id].size;
@@ -333,8 +331,7 @@
   }
 
   if (num_nonzeros_ > max_num_nonzeros_) {
-    std::unique_ptr<double[]> new_values =
-        std::make_unique<double[]>(num_nonzeros_);
+    auto new_values = std::make_unique<double[]>(num_nonzeros_);
     std::copy_n(values_.get(), old_num_nonzeros, new_values.get());
     values_ = std::move(new_values);
     max_num_nonzeros_ = num_nonzeros_;
diff --git a/internal/ceres/block_structure.cc b/internal/ceres/block_structure.cc
index f6767b4..cbc493b 100644
--- a/internal/ceres/block_structure.cc
+++ b/internal/ceres/block_structure.cc
@@ -30,6 +30,8 @@
 
 #include "ceres/block_structure.h"
 
+#include "glog/logging.h"
+
 namespace ceres::internal {
 
 bool CellLessThan(const Cell& lhs, const Cell& rhs) {
@@ -39,4 +41,28 @@
   return (lhs.block_id < rhs.block_id);
 }
 
+std::vector<Block> Tail(const std::vector<Block>& blocks, int n) {
+  CHECK_LE(n, blocks.size());
+  std::vector<Block> tail;
+  const int num_blocks = blocks.size();
+  const int start = num_blocks - n;
+
+  int position = 0;
+  tail.reserve(n);
+  for (int i = start; i < num_blocks; ++i) {
+    tail.emplace_back(blocks[i].size, position);
+    position += blocks[i].size;
+  }
+
+  return tail;
+}
+
+int SumSquaredSizes(const std::vector<Block>& blocks) {
+  int sum = 0;
+  for (const auto& b : blocks) {
+    sum += b.size * b.size;
+  }
+  return sum;
+}
+
 }  // namespace ceres::internal
diff --git a/internal/ceres/block_structure.h b/internal/ceres/block_structure.h
index 2ba6b32..9ee10af 100644
--- a/internal/ceres/block_structure.h
+++ b/internal/ceres/block_structure.h
@@ -49,15 +49,19 @@
 
 struct CERES_NO_EXPORT Block {
   Block() = default;
-  Block(int size_, int position_) : size(size_), position(position_) {}
+  Block(int size_, int position_) noexcept : size(size_), position(position_) {}
 
   BlockSize size{-1};
   int position{-1};  // Position along the row/column.
 };
 
+inline bool operator==(const Block& left, const Block& right) noexcept {
+  return (left.size == right.size) && (left.position == right.position);
+}
+
 struct CERES_NO_EXPORT Cell {
   Cell() = default;
-  Cell(int block_id_, int position_)
+  Cell(int block_id_, int position_) noexcept
       : block_id(block_id_), position(position_) {}
 
   // Column or row block id as the case maybe.
@@ -74,7 +78,7 @@
 
   // Construct a CompressedList with the cells containing num_cells
   // entries.
-  explicit CompressedList(int num_cells) : cells(num_cells) {}
+  explicit CompressedList(int num_cells) noexcept : cells(num_cells) {}
   Block block;
   std::vector<Cell> cells;
 };
@@ -168,6 +172,18 @@
   std::vector<CompressedColumn> cols;
 };
 
+inline int NumScalarEntries(const std::vector<Block>& blocks) {
+  if (blocks.empty()) {
+    return 0;
+  }
+
+  auto& block = blocks.back();
+  return block.position + block.size;
+}
+
+std::vector<Block> Tail(const std::vector<Block>& blocks, int n);
+int SumSquaredSizes(const std::vector<Block>& blocks);
+
 }  // namespace ceres::internal
 
 #endif  // CERES_INTERNAL_BLOCK_STRUCTURE_H_
diff --git a/internal/ceres/compressed_col_sparse_matrix_utils.cc b/internal/ceres/compressed_col_sparse_matrix_utils.cc
index cb6be55..8d94ec2 100644
--- a/internal/ceres/compressed_col_sparse_matrix_utils.cc
+++ b/internal/ceres/compressed_col_sparse_matrix_utils.cc
@@ -38,14 +38,13 @@
 
 namespace ceres::internal {
 
-using std::vector;
-
-void CompressedColumnScalarMatrixToBlockMatrix(const int* scalar_rows,
-                                               const int* scalar_cols,
-                                               const vector<int>& row_blocks,
-                                               const vector<int>& col_blocks,
-                                               vector<int>* block_rows,
-                                               vector<int>* block_cols) {
+void CompressedColumnScalarMatrixToBlockMatrix(
+    const int* scalar_rows,
+    const int* scalar_cols,
+    const std::vector<Block>& row_blocks,
+    const std::vector<Block>& col_blocks,
+    std::vector<int>* block_rows,
+    std::vector<int>* block_cols) {
   CHECK(block_rows != nullptr);
   CHECK(block_cols != nullptr);
   block_rows->clear();
@@ -53,12 +52,6 @@
   const int num_row_blocks = row_blocks.size();
   const int num_col_blocks = col_blocks.size();
 
-  vector<int> row_block_starts(num_row_blocks);
-  for (int i = 0, cursor = 0; i < num_row_blocks; ++i) {
-    row_block_starts[i] = cursor;
-    cursor += row_blocks[i];
-  }
-
   // This loop extracts the block sparsity of the scalar sparse matrix
   // It does so by iterating over the columns, but only considering
   // the columns corresponding to the first element of each column
@@ -70,48 +63,43 @@
   for (int col_block = 0; col_block < num_col_blocks; ++col_block) {
     int column_size = 0;
     for (int idx = scalar_cols[c]; idx < scalar_cols[c + 1]; ++idx) {
-      auto it = std::lower_bound(
-          row_block_starts.begin(), row_block_starts.end(), scalar_rows[idx]);
-      // Since we are using lower_bound, it will return the row id
-      // where the row block starts. For everything but the first row
-      // of the block, where these values will be the same, we can
-      // skip, as we only need the first row to detect the presence of
-      // the block.
+      auto it = std::lower_bound(row_blocks.begin(),
+                                 row_blocks.end(),
+                                 scalar_rows[idx],
+                                 [](const Block& block, double value) {
+                                   return block.position < value;
+                                 });
+      // Since we are using lower_bound, it will return the row id where the row
+      // block starts. For everything but the first row of the block, where
+      // these values will be the same, we can skip, as we only need the first
+      // row to detect the presence of the block.
       //
-      // For rows all but the first row in the last row block,
-      // lower_bound will return row_block_starts.end(), but those can
-      // be skipped like the rows in other row blocks too.
-      if (it == row_block_starts.end() || *it != scalar_rows[idx]) {
+      // For rows all but the first row in the last row block, lower_bound will
+      // return row_blocks_.end(), but those can be skipped like the rows in
+      // other row blocks too.
+      if (it == row_blocks.end() || it->position != scalar_rows[idx]) {
         continue;
       }
 
-      block_rows->push_back(it - row_block_starts.begin());
+      block_rows->push_back(it - row_blocks.begin());
       ++column_size;
     }
     block_cols->push_back(block_cols->back() + column_size);
-    c += col_blocks[col_block];
+    c += col_blocks[col_block].size;
   }
 }
 
-void BlockOrderingToScalarOrdering(const vector<int>& blocks,
-                                   const vector<int>& block_ordering,
-                                   vector<int>* scalar_ordering) {
+void BlockOrderingToScalarOrdering(const std::vector<Block>& blocks,
+                                   const std::vector<int>& block_ordering,
+                                   std::vector<int>* scalar_ordering) {
   CHECK_EQ(blocks.size(), block_ordering.size());
   const int num_blocks = blocks.size();
-
-  // block_starts = [0, block1, block1 + block2 ..]
-  vector<int> block_starts(num_blocks);
-  for (int i = 0, cursor = 0; i < num_blocks; ++i) {
-    block_starts[i] = cursor;
-    cursor += blocks[i];
-  }
-
-  scalar_ordering->resize(block_starts.back() + blocks.back());
+  scalar_ordering->resize(NumScalarEntries(blocks));
   int cursor = 0;
   for (int i = 0; i < num_blocks; ++i) {
     const int block_id = block_ordering[i];
-    const int block_size = blocks[block_id];
-    int block_position = block_starts[block_id];
+    const int block_size = blocks[block_id].size;
+    int block_position = blocks[block_id].position;
     for (int j = 0; j < block_size; ++j) {
       (*scalar_ordering)[cursor++] = block_position++;
     }
diff --git a/internal/ceres/compressed_col_sparse_matrix_utils.h b/internal/ceres/compressed_col_sparse_matrix_utils.h
index 608150e..fbb0abf 100644
--- a/internal/ceres/compressed_col_sparse_matrix_utils.h
+++ b/internal/ceres/compressed_col_sparse_matrix_utils.h
@@ -34,6 +34,7 @@
 #include <algorithm>
 #include <vector>
 
+#include "ceres/block_structure.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
 
@@ -52,8 +53,8 @@
 CERES_NO_EXPORT void CompressedColumnScalarMatrixToBlockMatrix(
     const int* scalar_rows,
     const int* scalar_cols,
-    const std::vector<int>& row_blocks,
-    const std::vector<int>& col_blocks,
+    const std::vector<Block>& row_blocks,
+    const std::vector<Block>& col_blocks,
     std::vector<int>* block_rows,
     std::vector<int>* block_cols);
 
@@ -61,7 +62,7 @@
 // the corresponding "scalar" ordering, where the scalar ordering of
 // size sum(blocks).
 CERES_NO_EXPORT void BlockOrderingToScalarOrdering(
-    const std::vector<int>& blocks,
+    const std::vector<Block>& blocks,
     const std::vector<int>& block_ordering,
     std::vector<int>* scalar_ordering);
 
diff --git a/internal/ceres/compressed_col_sparse_matrix_utils_test.cc b/internal/ceres/compressed_col_sparse_matrix_utils_test.cc
index 3b7f9de..ef42ee2 100644
--- a/internal/ceres/compressed_col_sparse_matrix_utils_test.cc
+++ b/internal/ceres/compressed_col_sparse_matrix_utils_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -39,46 +39,22 @@
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
-namespace ceres {
-namespace internal {
-
-using std::vector;
+namespace ceres::internal {
 
 TEST(_, BlockPermutationToScalarPermutation) {
-  vector<int> blocks;
   //  Block structure
   //  0  --1-  ---2---  ---3---  4
   // [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
-  blocks.push_back(1);
-  blocks.push_back(2);
-  blocks.push_back(3);
-  blocks.push_back(3);
-  blocks.push_back(1);
-
+  std::vector<Block> blocks{{1, 0}, {2, 1}, {3, 3}, {3, 6}, {1, 9}};
   // Block ordering
   // [1, 0, 2, 4, 5]
-  vector<int> block_ordering;
-  block_ordering.push_back(1);
-  block_ordering.push_back(0);
-  block_ordering.push_back(2);
-  block_ordering.push_back(4);
-  block_ordering.push_back(3);
+  std::vector<int> block_ordering{{1, 0, 2, 4, 3}};
 
   // Expected ordering
   // [1, 2, 0, 3, 4, 5, 9, 6, 7, 8]
-  vector<int> expected_scalar_ordering;
-  expected_scalar_ordering.push_back(1);
-  expected_scalar_ordering.push_back(2);
-  expected_scalar_ordering.push_back(0);
-  expected_scalar_ordering.push_back(3);
-  expected_scalar_ordering.push_back(4);
-  expected_scalar_ordering.push_back(5);
-  expected_scalar_ordering.push_back(9);
-  expected_scalar_ordering.push_back(6);
-  expected_scalar_ordering.push_back(7);
-  expected_scalar_ordering.push_back(8);
+  std::vector<int> expected_scalar_ordering{{1, 2, 0, 3, 4, 5, 9, 6, 7, 8}};
 
-  vector<int> scalar_ordering;
+  std::vector<int> scalar_ordering;
   BlockOrderingToScalarOrdering(blocks, block_ordering, &scalar_ordering);
   EXPECT_EQ(scalar_ordering.size(), expected_scalar_ordering.size());
   for (int i = 0; i < expected_scalar_ordering.size(); ++i) {
@@ -86,19 +62,17 @@
   }
 }
 
-static void FillBlock(const vector<int>& row_blocks,
-                      const vector<int>& col_blocks,
+static void FillBlock(const std::vector<Block>& row_blocks,
+                      const std::vector<Block>& col_blocks,
                       const int row_block_id,
                       const int col_block_id,
-                      vector<Eigen::Triplet<double>>* triplets) {
-  const int row_offset =
-      std::accumulate(&row_blocks[0], &row_blocks[row_block_id], 0);
-  const int col_offset =
-      std::accumulate(&col_blocks[0], &col_blocks[col_block_id], 0);
-  for (int r = 0; r < row_blocks[row_block_id]; ++r) {
-    for (int c = 0; c < col_blocks[col_block_id]; ++c) {
+                      std::vector<Eigen::Triplet<double>>* triplets) {
+  for (int r = 0; r < row_blocks[row_block_id].size; ++r) {
+    for (int c = 0; c < col_blocks[col_block_id].size; ++c) {
       triplets->push_back(
-          Eigen::Triplet<double>(row_offset + r, col_offset + c, 1.0));
+          Eigen::Triplet<double>(row_blocks[row_block_id].position + r,
+                                 col_blocks[col_block_id].position + c,
+                                 1.0));
     }
   }
 }
@@ -112,23 +86,13 @@
   // [2]  x x
   // num_nonzeros = 1 + 3 + 4 + 4 + 1 + 2 = 15
 
-  vector<int> col_blocks;
-  col_blocks.push_back(1);
-  col_blocks.push_back(2);
-  col_blocks.push_back(3);
-  col_blocks.push_back(2);
+  std::vector<Block> col_blocks{{1, 0}, {2, 1}, {3, 3}, {2, 5}};
+  const int num_cols = NumScalarEntries(col_blocks);
 
-  vector<int> row_blocks;
-  row_blocks.push_back(1);
-  row_blocks.push_back(2);
-  row_blocks.push_back(2);
+  std::vector<Block> row_blocks{{1, 0}, {2, 1}, {2, 3}};
+  const int num_rows = NumScalarEntries(row_blocks);
 
-  const int num_rows =
-      std::accumulate(row_blocks.begin(), row_blocks.end(), 0.0);
-  const int num_cols =
-      std::accumulate(col_blocks.begin(), col_blocks.end(), 0.0);
-
-  vector<Eigen::Triplet<double>> triplets;
+  std::vector<Eigen::Triplet<double>> triplets;
   FillBlock(row_blocks, col_blocks, 0, 0, &triplets);
   FillBlock(row_blocks, col_blocks, 2, 0, &triplets);
   FillBlock(row_blocks, col_blocks, 1, 1, &triplets);
@@ -138,23 +102,11 @@
   Eigen::SparseMatrix<double> sparse_matrix(num_rows, num_cols);
   sparse_matrix.setFromTriplets(triplets.begin(), triplets.end());
 
-  vector<int> expected_compressed_block_rows;
-  expected_compressed_block_rows.push_back(0);
-  expected_compressed_block_rows.push_back(2);
-  expected_compressed_block_rows.push_back(1);
-  expected_compressed_block_rows.push_back(2);
-  expected_compressed_block_rows.push_back(0);
-  expected_compressed_block_rows.push_back(1);
+  const std::vector<int> expected_compressed_block_rows{{0, 2, 1, 2, 0, 1}};
+  const std::vector<int> expected_compressed_block_cols{{0, 2, 4, 5, 6}};
 
-  vector<int> expected_compressed_block_cols;
-  expected_compressed_block_cols.push_back(0);
-  expected_compressed_block_cols.push_back(2);
-  expected_compressed_block_cols.push_back(4);
-  expected_compressed_block_cols.push_back(5);
-  expected_compressed_block_cols.push_back(6);
-
-  vector<int> compressed_block_rows;
-  vector<int> compressed_block_cols;
+  std::vector<int> compressed_block_rows;
+  std::vector<int> compressed_block_cols;
   CompressedColumnScalarMatrixToBlockMatrix(sparse_matrix.innerIndexPtr(),
                                             sparse_matrix.outerIndexPtr(),
                                             row_blocks,
@@ -198,9 +150,9 @@
     cols[4] = 7;
   }
 
-  vector<int> cols;
-  vector<int> rows;
-  vector<double> values;
+  std::vector<int> cols;
+  std::vector<int> rows;
+  std::vector<double> values;
 };
 
 TEST_F(SolveUpperTriangularTest, SolveInPlace) {
@@ -245,5 +197,4 @@
   }
 }
 
-}  // namespace internal
-}  // namespace ceres
+}  // namespace ceres::internal
diff --git a/internal/ceres/compressed_row_jacobian_writer.cc b/internal/ceres/compressed_row_jacobian_writer.cc
index c409872..7112554 100644
--- a/internal/ceres/compressed_row_jacobian_writer.cc
+++ b/internal/ceres/compressed_row_jacobian_writer.cc
@@ -45,42 +45,41 @@
 #include "ceres/scratch_evaluate_preparer.h"
 
 namespace ceres::internal {
-
-using std::adjacent_find;
-using std::make_pair;
-using std::pair;
-using std::vector;
-
 void CompressedRowJacobianWriter::PopulateJacobianRowAndColumnBlockVectors(
     const Program* program, CompressedRowSparseMatrix* jacobian) {
-  const vector<ParameterBlock*>& parameter_blocks = program->parameter_blocks();
-  vector<int>& col_blocks = *(jacobian->mutable_col_blocks());
+  const auto& parameter_blocks = program->parameter_blocks();
+  auto& col_blocks = *(jacobian->mutable_col_blocks());
   col_blocks.resize(parameter_blocks.size());
+  int col_pos = 0;
   for (int i = 0; i < parameter_blocks.size(); ++i) {
-    col_blocks[i] = parameter_blocks[i]->TangentSize();
+    col_blocks[i].size = parameter_blocks[i]->TangentSize();
+    col_blocks[i].position = col_pos;
+    col_pos += col_blocks[i].size;
   }
 
-  const vector<ResidualBlock*>& residual_blocks = program->residual_blocks();
-  vector<int>& row_blocks = *(jacobian->mutable_row_blocks());
+  const auto& residual_blocks = program->residual_blocks();
+  auto& row_blocks = *(jacobian->mutable_row_blocks());
   row_blocks.resize(residual_blocks.size());
+  int row_pos = 0;
   for (int i = 0; i < residual_blocks.size(); ++i) {
-    row_blocks[i] = residual_blocks[i]->NumResiduals();
+    row_blocks[i].size = residual_blocks[i]->NumResiduals();
+    row_blocks[i].position = row_pos;
+    row_pos += row_blocks[i].size;
   }
 }
 
 void CompressedRowJacobianWriter::GetOrderedParameterBlocks(
     const Program* program,
     int residual_id,
-    vector<pair<int, int>>* evaluated_jacobian_blocks) {
-  const ResidualBlock* residual_block = program->residual_blocks()[residual_id];
+    std::vector<std::pair<int, int>>* evaluated_jacobian_blocks) {
+  auto residual_block = program->residual_blocks()[residual_id];
   const int num_parameter_blocks = residual_block->NumParameterBlocks();
 
   for (int j = 0; j < num_parameter_blocks; ++j) {
-    const ParameterBlock* parameter_block =
-        residual_block->parameter_blocks()[j];
+    auto parameter_block = residual_block->parameter_blocks()[j];
     if (!parameter_block->IsConstant()) {
       evaluated_jacobian_blocks->push_back(
-          make_pair(parameter_block->index(), j));
+          std::make_pair(parameter_block->index(), j));
     }
   }
   std::sort(evaluated_jacobian_blocks->begin(),
@@ -89,7 +88,7 @@
 
 std::unique_ptr<SparseMatrix> CompressedRowJacobianWriter::CreateJacobian()
     const {
-  const vector<ResidualBlock*>& residual_blocks = program_->residual_blocks();
+  const auto& residual_blocks = program_->residual_blocks();
 
   int total_num_residuals = program_->NumResiduals();
   int total_num_effective_parameters = program_->NumEffectiveParameters();
@@ -100,7 +99,7 @@
     const int num_residuals = residual_block->NumResiduals();
     const int num_parameter_blocks = residual_block->NumParameterBlocks();
     for (int j = 0; j < num_parameter_blocks; ++j) {
-      ParameterBlock* parameter_block = residual_block->parameter_blocks()[j];
+      auto parameter_block = residual_block->parameter_blocks()[j];
       if (!parameter_block->IsConstant()) {
         num_jacobian_nonzeros += num_residuals * parameter_block->TangentSize();
       }
@@ -111,14 +110,14 @@
   // Allocate more space than needed to store the jacobian so that when the LM
   // algorithm adds the diagonal, no reallocation is necessary. This reduces
   // peak memory usage significantly.
-  std::unique_ptr<CompressedRowSparseMatrix> jacobian =
-      std::make_unique<CompressedRowSparseMatrix>(
-          total_num_residuals,
-          total_num_effective_parameters,
-          num_jacobian_nonzeros + total_num_effective_parameters);
+  auto jacobian = std::make_unique<CompressedRowSparseMatrix>(
+      total_num_residuals,
+      total_num_effective_parameters,
+      num_jacobian_nonzeros + total_num_effective_parameters);
 
-  // At this stage, the CompressedRowSparseMatrix is an invalid state. But this
-  // seems to be the only way to construct it without doing a memory copy.
+  // At this stage, the CompressedRowSparseMatrix is an invalid state. But
+  // this seems to be the only way to construct it without doing a memory
+  // copy.
   int* rows = jacobian->mutable_rows();
   int* cols = jacobian->mutable_cols();
 
@@ -130,9 +129,9 @@
     // Count the number of derivatives for a row of this residual block and
     // build a list of active parameter block indices.
     int num_derivatives = 0;
-    vector<int> parameter_indices;
+    std::vector<int> parameter_indices;
     for (int j = 0; j < num_parameter_blocks; ++j) {
-      ParameterBlock* parameter_block = residual_block->parameter_blocks()[j];
+      auto parameter_block = residual_block->parameter_blocks()[j];
       if (!parameter_block->IsConstant()) {
         parameter_indices.push_back(parameter_block->index());
         num_derivatives += parameter_block->TangentSize();
@@ -140,12 +139,12 @@
     }
 
     // Sort the parameters by their position in the state vector.
-    sort(parameter_indices.begin(), parameter_indices.end());
+    std::sort(parameter_indices.begin(), parameter_indices.end());
     if (adjacent_find(parameter_indices.begin(), parameter_indices.end()) !=
         parameter_indices.end()) {
       std::string parameter_block_description;
       for (int j = 0; j < num_parameter_blocks; ++j) {
-        ParameterBlock* parameter_block = residual_block->parameter_blocks()[j];
+        auto parameter_block = residual_block->parameter_blocks()[j];
         parameter_block_description += parameter_block->ToString() + "\n";
       }
       LOG(FATAL) << "Ceres internal error: "
@@ -167,15 +166,13 @@
     // values are updated.
     int col_pos = 0;
     for (int parameter_index : parameter_indices) {
-      ParameterBlock* parameter_block =
-          program_->parameter_blocks()[parameter_index];
+      auto parameter_block = program_->parameter_blocks()[parameter_index];
       const int parameter_block_size = parameter_block->TangentSize();
 
       for (int r = 0; r < num_residuals; ++r) {
         // This is the position in the values array of the jacobian where this
         // row of the jacobian block should go.
         const int column_block_begin = rows[row_pos + r] + col_pos;
-
         for (int c = 0; c < parameter_block_size; ++c) {
           cols[column_block_begin + c] = parameter_block->delta_offset() + c;
         }
@@ -200,11 +197,10 @@
   double* jacobian_values = jacobian->mutable_values();
   const int* jacobian_rows = jacobian->rows();
 
-  const ResidualBlock* residual_block =
-      program_->residual_blocks()[residual_id];
+  auto residual_block = program_->residual_blocks()[residual_id];
   const int num_residuals = residual_block->NumResiduals();
 
-  vector<pair<int, int>> evaluated_jacobian_blocks;
+  std::vector<std::pair<int, int>> evaluated_jacobian_blocks;
   GetOrderedParameterBlocks(program_, residual_id, &evaluated_jacobian_blocks);
 
   // Where in the current row does the jacobian for a parameter block begin.
@@ -213,7 +209,7 @@
   // Iterate over the jacobian blocks in increasing order of their
   // positions in the reduced parameter vector.
   for (auto& evaluated_jacobian_block : evaluated_jacobian_blocks) {
-    const ParameterBlock* parameter_block =
+    auto parameter_block =
         program_->parameter_blocks()[evaluated_jacobian_block.first];
     const int argument = evaluated_jacobian_block.second;
     const int parameter_block_size = parameter_block->TangentSize();
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index f2c7747..057adfd 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -43,9 +43,6 @@
 #include "glog/logging.h"
 
 namespace ceres::internal {
-
-using std::vector;
-
 namespace {
 
 // Helper functor used by the constructor for reordering the contents
@@ -206,7 +203,7 @@
   }
 
   // index is the list of indices into the TripletSparseMatrix input.
-  vector<int> index(input.num_nonzeros(), 0);
+  std::vector<int> index(input.num_nonzeros(), 0);
   for (int i = 0; i < input.num_nonzeros(); ++i) {
     index[i] = i;
   }
@@ -221,9 +218,8 @@
               input.num_nonzeros() * sizeof(int) +     // NOLINT
               input.num_nonzeros() * sizeof(double));  // NOLINT
 
-  std::unique_ptr<CompressedRowSparseMatrix> output =
-      std::make_unique<CompressedRowSparseMatrix>(
-          num_rows, num_cols, input.num_nonzeros());
+  auto output = std::make_unique<CompressedRowSparseMatrix>(
+      num_rows, num_cols, input.num_nonzeros());
 
   if (num_rows == 0) {
     // No data to copy.
@@ -457,7 +453,7 @@
   int num_row_blocks = 0;
   int num_rows = 0;
   while (num_row_blocks < row_blocks_.size() && num_rows < num_rows_) {
-    num_rows += row_blocks_[num_row_blocks];
+    num_rows += row_blocks_[num_row_blocks].size;
     ++num_row_blocks;
   }
 
@@ -545,17 +541,15 @@
 
 std::unique_ptr<CompressedRowSparseMatrix>
 CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
-    const double* diagonal, const vector<int>& blocks) {
-  int num_rows = 0;
+    const double* diagonal, const std::vector<Block>& blocks) {
+  const int num_rows = NumScalarEntries(blocks);
   int num_nonzeros = 0;
-  for (int block_size : blocks) {
-    num_rows += block_size;
-    num_nonzeros += block_size * block_size;
+  for (auto& block : blocks) {
+    num_nonzeros += block.size * block.size;
   }
 
-  std::unique_ptr<CompressedRowSparseMatrix> matrix =
-      std::make_unique<CompressedRowSparseMatrix>(
-          num_rows, num_rows, num_nonzeros);
+  auto matrix = std::make_unique<CompressedRowSparseMatrix>(
+      num_rows, num_rows, num_nonzeros);
 
   int* rows = matrix->mutable_rows();
   int* cols = matrix->mutable_cols();
@@ -564,15 +558,15 @@
 
   int idx_cursor = 0;
   int col_cursor = 0;
-  for (int block_size : blocks) {
-    for (int r = 0; r < block_size; ++r) {
+  for (auto& block : blocks) {
+    for (int r = 0; r < block.size; ++r) {
       *(rows++) = idx_cursor;
       values[idx_cursor + r] = diagonal[col_cursor + r];
-      for (int c = 0; c < block_size; ++c, ++idx_cursor) {
+      for (int c = 0; c < block.size; ++c, ++idx_cursor) {
         *(cols++) = col_cursor + c;
       }
     }
-    col_cursor += block_size;
+    col_cursor += block.size;
   }
   *rows = idx_cursor;
 
@@ -586,9 +580,8 @@
 
 std::unique_ptr<CompressedRowSparseMatrix>
 CompressedRowSparseMatrix::Transpose() const {
-  std::unique_ptr<CompressedRowSparseMatrix> transpose =
-      std::make_unique<CompressedRowSparseMatrix>(
-          num_cols_, num_rows_, num_nonzeros());
+  auto transpose = std::make_unique<CompressedRowSparseMatrix>(
+      num_cols_, num_rows_, num_nonzeros());
 
   switch (storage_type_) {
     case StorageType::UNSYMMETRIC:
@@ -649,9 +642,9 @@
   CHECK_GT(options.block_density, 0.0);
   CHECK_LE(options.block_density, 1.0);
 
-  vector<int> row_blocks;
+  std::vector<Block> row_blocks;
   row_blocks.reserve(options.num_row_blocks);
-  vector<int> col_blocks;
+  std::vector<Block> col_blocks;
   col_blocks.reserve(options.num_col_blocks);
 
   std::uniform_int_distribution<int> col_distribution(
@@ -662,25 +655,29 @@
   std::normal_distribution<double> standard_normal_distribution;
 
   // Generate the row block structure.
+  int row_pos = 0;
   for (int i = 0; i < options.num_row_blocks; ++i) {
     // Generate a random integer in [min_row_block_size, max_row_block_size]
-    row_blocks.push_back(row_distribution(prng));
+    row_blocks.emplace_back(row_distribution(prng), row_pos);
+    row_pos += row_blocks.back().size;
   }
 
   if (options.storage_type == StorageType::UNSYMMETRIC) {
     // Generate the col block structure.
+    int col_pos = 0;
     for (int i = 0; i < options.num_col_blocks; ++i) {
       // Generate a random integer in [min_col_block_size, max_col_block_size]
-      col_blocks.push_back(col_distribution(prng));
+      col_blocks.emplace_back(col_distribution(prng), col_pos);
+      col_pos += col_blocks.back().size;
     }
   } else {
     // Symmetric matrices (LOWER_TRIANGULAR or UPPER_TRIANGULAR);
     col_blocks = row_blocks;
   }
 
-  vector<int> tsm_rows;
-  vector<int> tsm_cols;
-  vector<double> tsm_values;
+  std::vector<int> tsm_rows;
+  std::vector<int> tsm_cols;
+  std::vector<double> tsm_values;
 
   // For ease of construction, we are going to generate the
   // CompressedRowSparseMatrix by generating it as a
@@ -703,7 +700,7 @@
              (r > c)) ||
             ((options.storage_type == StorageType::LOWER_TRIANGULAR) &&
              (r < c))) {
-          col_block_begin += col_blocks[c];
+          col_block_begin += col_blocks[c].size;
           continue;
         }
 
@@ -715,8 +712,8 @@
           // If the matrix is symmetric, then we take care to generate
           // symmetric diagonal blocks.
           if (options.storage_type == StorageType::UNSYMMETRIC || r != c) {
-            AddRandomBlock(row_blocks[r],
-                           col_blocks[c],
+            AddRandomBlock(row_blocks[r].size,
+                           col_blocks[c].size,
                            row_block_begin,
                            col_block_begin,
                            randn,
@@ -724,7 +721,7 @@
                            &tsm_cols,
                            &tsm_values);
           } else {
-            AddSymmetricRandomBlock(row_blocks[r],
+            AddSymmetricRandomBlock(row_blocks[r].size,
                                     row_block_begin,
                                     randn,
                                     &tsm_rows,
@@ -732,20 +729,18 @@
                                     &tsm_values);
           }
         }
-        col_block_begin += col_blocks[c];
+        col_block_begin += col_blocks[c].size;
       }
-      row_block_begin += row_blocks[r];
+      row_block_begin += row_blocks[r].size;
     }
   }
 
-  const int num_rows = std::accumulate(row_blocks.begin(), row_blocks.end(), 0);
-  const int num_cols = std::accumulate(col_blocks.begin(), col_blocks.end(), 0);
+  const int num_rows = NumScalarEntries(row_blocks);
+  const int num_cols = NumScalarEntries(col_blocks);
   const bool kDoNotTranspose = false;
-  std::unique_ptr<CompressedRowSparseMatrix> matrix =
-      CompressedRowSparseMatrix::FromTripletSparseMatrix(
-          TripletSparseMatrix(
-              num_rows, num_cols, tsm_rows, tsm_cols, tsm_values),
-          kDoNotTranspose);
+  auto matrix = CompressedRowSparseMatrix::FromTripletSparseMatrix(
+      TripletSparseMatrix(num_rows, num_cols, tsm_rows, tsm_cols, tsm_values),
+      kDoNotTranspose);
   (*matrix->mutable_row_blocks()) = row_blocks;
   (*matrix->mutable_col_blocks()) = col_blocks;
   matrix->set_storage_type(options.storage_type);
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index d98f61f..8668d06 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -35,6 +35,7 @@
 #include <random>
 #include <vector>
 
+#include "ceres/block_structure.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
 #include "ceres/sparse_matrix.h"
@@ -144,17 +145,17 @@
     storage_type_ = storage_type;
   }
 
-  const std::vector<int>& row_blocks() const { return row_blocks_; }
-  std::vector<int>* mutable_row_blocks() { return &row_blocks_; }
+  const std::vector<Block>& row_blocks() const { return row_blocks_; }
+  std::vector<Block>* mutable_row_blocks() { return &row_blocks_; }
 
-  const std::vector<int>& col_blocks() const { return col_blocks_; }
-  std::vector<int>* mutable_col_blocks() { return &col_blocks_; }
+  const std::vector<Block>& col_blocks() const { return col_blocks_; }
+  std::vector<Block>* mutable_col_blocks() { return &col_blocks_; }
 
   // Create a block diagonal CompressedRowSparseMatrix with the given
   // block structure. The individual blocks are assumed to be laid out
   // contiguously in the diagonal array, one block at a time.
   static std::unique_ptr<CompressedRowSparseMatrix> CreateBlockDiagonalMatrix(
-      const double* diagonal, const std::vector<int>& blocks);
+      const double* diagonal, const std::vector<Block>& blocks);
 
   // Options struct to control the generation of random block sparse
   // matrices in compressed row sparse format.
@@ -214,8 +215,8 @@
   // optional information for use by algorithms operating on the
   // matrix. The class itself does not make use of this information in
   // any way.
-  std::vector<int> row_blocks_;
-  std::vector<int> col_blocks_;
+  std::vector<Block> row_blocks_;
+  std::vector<Block> col_blocks_;
 };
 
 inline std::ostream& operator<<(std::ostream& s,
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index 718a659..82c982b 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -35,6 +35,7 @@
 #include <numeric>
 #include <random>
 #include <string>
+#include <vector>
 
 #include "Eigen/SparseCore"
 #include "ceres/casts.h"
@@ -45,10 +46,7 @@
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
-namespace ceres {
-namespace internal {
-
-using std::vector;
+namespace ceres::internal {
 
 static void CompareMatrices(const SparseMatrix* a, const SparseMatrix* b) {
   EXPECT_EQ(a->num_rows(), b->num_rows());
@@ -73,8 +71,7 @@
 class CompressedRowSparseMatrixTest : public ::testing::Test {
  protected:
   void SetUp() final {
-    std::unique_ptr<LinearLeastSquaresProblem> problem =
-        CreateLinearLeastSquaresProblemFromId(1);
+    auto problem = CreateLinearLeastSquaresProblemFromId(1);
 
     CHECK(problem != nullptr);
 
@@ -84,13 +81,16 @@
     num_rows = tsm->num_rows();
     num_cols = tsm->num_cols();
 
-    vector<int>* row_blocks = crsm->mutable_row_blocks();
+    std::vector<Block>* row_blocks = crsm->mutable_row_blocks();
     row_blocks->resize(num_rows);
-    std::fill(row_blocks->begin(), row_blocks->end(), 1);
-
-    vector<int>* col_blocks = crsm->mutable_col_blocks();
+    for (int i = 0; i < row_blocks->size(); ++i) {
+      (*row_blocks)[i] = Block(1, i);
+    }
+    std::vector<Block>* col_blocks = crsm->mutable_col_blocks();
     col_blocks->resize(num_cols);
-    std::fill(col_blocks->begin(), col_blocks->end(), 1);
+    for (int i = 0; i < col_blocks->size(); ++i) {
+      (*col_blocks)[i] = Block(1, i);
+    }
   }
 
   int num_rows;
@@ -133,7 +133,7 @@
     tsm_appendage.Resize(i, num_cols);
 
     tsm->AppendRows(tsm_appendage);
-    std::unique_ptr<CompressedRowSparseMatrix> crsm_appendage =
+    auto crsm_appendage =
         CompressedRowSparseMatrix::FromTripletSparseMatrix(tsm_appendage);
 
     crsm->AppendRows(*crsm_appendage);
@@ -144,35 +144,33 @@
 TEST_F(CompressedRowSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) {
   int num_diagonal_rows = crsm->num_cols();
 
-  std::unique_ptr<double[]> diagonal =
-      std::make_unique<double[]>(num_diagonal_rows);
+  auto diagonal = std::make_unique<double[]>(num_diagonal_rows);
   for (int i = 0; i < num_diagonal_rows; ++i) {
     diagonal[i] = i;
   }
 
-  vector<int> row_and_column_blocks;
-  row_and_column_blocks.push_back(1);
-  row_and_column_blocks.push_back(2);
-  row_and_column_blocks.push_back(2);
+  std::vector<Block> row_and_column_blocks;
+  row_and_column_blocks.emplace_back(1, 0);
+  row_and_column_blocks.emplace_back(2, 1);
+  row_and_column_blocks.emplace_back(2, 3);
 
-  const vector<int> pre_row_blocks = crsm->row_blocks();
-  const vector<int> pre_col_blocks = crsm->col_blocks();
+  const std::vector<Block> pre_row_blocks = crsm->row_blocks();
+  const std::vector<Block> pre_col_blocks = crsm->col_blocks();
 
-  std::unique_ptr<CompressedRowSparseMatrix> appendage =
-      CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
-          diagonal.get(), row_and_column_blocks);
+  auto appendage = CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
+      diagonal.get(), row_and_column_blocks);
 
   crsm->AppendRows(*appendage);
 
-  const vector<int> post_row_blocks = crsm->row_blocks();
-  const vector<int> post_col_blocks = crsm->col_blocks();
+  const std::vector<Block> post_row_blocks = crsm->row_blocks();
+  const std::vector<Block> post_col_blocks = crsm->col_blocks();
 
-  vector<int> expected_row_blocks = pre_row_blocks;
+  std::vector<Block> expected_row_blocks = pre_row_blocks;
   expected_row_blocks.insert(expected_row_blocks.end(),
                              row_and_column_blocks.begin(),
                              row_and_column_blocks.end());
 
-  vector<int> expected_col_blocks = pre_col_blocks;
+  std::vector<Block> expected_col_blocks = pre_col_blocks;
 
   EXPECT_EQ(expected_row_blocks, crsm->row_blocks());
   EXPECT_EQ(expected_col_blocks, crsm->col_blocks());
@@ -212,19 +210,18 @@
 }
 
 TEST(CompressedRowSparseMatrix, CreateBlockDiagonalMatrix) {
-  vector<int> blocks;
-  blocks.push_back(1);
-  blocks.push_back(2);
-  blocks.push_back(2);
+  std::vector<Block> blocks;
+  blocks.emplace_back(1, 0);
+  blocks.emplace_back(2, 1);
+  blocks.emplace_back(2, 3);
 
   Vector diagonal(5);
   for (int i = 0; i < 5; ++i) {
     diagonal(i) = i + 1;
   }
 
-  std::unique_ptr<CompressedRowSparseMatrix> matrix =
-      CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(diagonal.data(),
-                                                           blocks);
+  auto matrix = CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
+      diagonal.data(), blocks);
 
   EXPECT_EQ(matrix->num_rows(), 5);
   EXPECT_EQ(matrix->num_cols(), 5);
@@ -272,10 +269,10 @@
   int* rows = matrix.mutable_rows();
   int* cols = matrix.mutable_cols();
   double* values = matrix.mutable_values();
-  matrix.mutable_row_blocks()->push_back(3);
-  matrix.mutable_row_blocks()->push_back(3);
-  matrix.mutable_col_blocks()->push_back(4);
-  matrix.mutable_col_blocks()->push_back(2);
+  matrix.mutable_row_blocks()->emplace_back(3, 0);
+  matrix.mutable_row_blocks()->emplace_back(3, 3);
+  matrix.mutable_col_blocks()->emplace_back(4, 0);
+  matrix.mutable_col_blocks()->emplace_back(2, 4);
 
   rows[0] = 0;
   cols[0] = 1;
@@ -307,7 +304,7 @@
 
   std::copy(values, values + 17, cols);
 
-  std::unique_ptr<CompressedRowSparseMatrix> transpose = matrix.Transpose();
+  auto transpose = matrix.Transpose();
 
   ASSERT_EQ(transpose->row_blocks().size(), matrix.col_blocks().size());
   for (int i = 0; i < transpose->row_blocks().size(); ++i) {
@@ -336,10 +333,8 @@
 
   const int kNumTrials = 10;
   for (int i = 0; i < kNumTrials; ++i) {
-    std::unique_ptr<TripletSparseMatrix> tsm =
-        TripletSparseMatrix::CreateRandomMatrix(options, prng);
-    std::unique_ptr<CompressedRowSparseMatrix> crsm =
-        CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm);
+    auto tsm = TripletSparseMatrix::CreateRandomMatrix(options, prng);
+    auto crsm = CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm);
 
     Matrix expected;
     tsm->ToDenseMatrix(&expected);
@@ -363,9 +358,8 @@
 
   const int kNumTrials = 10;
   for (int i = 0; i < kNumTrials; ++i) {
-    std::unique_ptr<TripletSparseMatrix> tsm =
-        TripletSparseMatrix::CreateRandomMatrix(options, prng);
-    std::unique_ptr<CompressedRowSparseMatrix> crsm =
+    auto tsm = TripletSparseMatrix::CreateRandomMatrix(options, prng);
+    auto crsm =
         CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm);
 
     Matrix tmp;
@@ -422,7 +416,7 @@
       options.max_row_block_size = kMaxBlockSize;
       options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
-      std::unique_ptr<CompressedRowSparseMatrix> matrix =
+      auto matrix =
           CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
       const int num_rows = matrix->num_rows();
       const int num_cols = matrix->num_cols();
@@ -492,7 +486,7 @@
       options.max_row_block_size = kMaxBlockSize;
       options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
-      std::unique_ptr<CompressedRowSparseMatrix> matrix =
+      auto matrix =
           CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
       const int num_rows = matrix->num_rows();
       const int num_cols = matrix->num_cols();
@@ -562,7 +556,7 @@
       options.max_row_block_size = kMaxBlockSize;
       options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
-      std::unique_ptr<CompressedRowSparseMatrix> matrix =
+      auto matrix =
           CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
       const int num_cols = matrix->num_cols();
 
@@ -607,5 +601,4 @@
 
 // TODO(sameeragarwal) Add tests for the random matrix creation methods.
 
-}  // namespace internal
-}  // namespace ceres
+}  // namespace ceres::internal
diff --git a/internal/ceres/implicit_schur_complement_test.cc b/internal/ceres/implicit_schur_complement_test.cc
index f7abca1..f0bbf6d 100644
--- a/internal/ceres/implicit_schur_complement_test.cc
+++ b/internal/ceres/implicit_schur_complement_test.cc
@@ -56,8 +56,7 @@
 class ImplicitSchurComplementTest : public ::testing::Test {
  protected:
   void SetUp() final {
-    std::unique_ptr<LinearLeastSquaresProblem> problem =
-        CreateLinearLeastSquaresProblemFromId(2);
+    auto problem = CreateLinearLeastSquaresProblemFromId(2);
 
     CHECK(problem != nullptr);
     A_.reset(down_cast<BlockSparseMatrix*>(problem->A.release()));
@@ -75,11 +74,7 @@
                                       Vector* solution) {
     const CompressedRowBlockStructure* bs = A_->block_structure();
     const int num_col_blocks = bs->cols.size();
-    std::vector<int> blocks(num_col_blocks - num_eliminate_blocks_, 0);
-    for (int i = num_eliminate_blocks_; i < num_col_blocks; ++i) {
-      blocks[i - num_eliminate_blocks_] = bs->cols[i].size;
-    }
-
+    auto blocks = Tail(bs->cols, num_col_blocks - num_eliminate_blocks_);
     BlockRandomAccessDenseMatrix blhs(blocks);
     const int num_schur_rows = blhs.num_rows();
 
diff --git a/internal/ceres/inner_product_computer.cc b/internal/ceres/inner_product_computer.cc
index f2d5457..a1d6832 100644
--- a/internal/ceres/inner_product_computer.cc
+++ b/internal/ceres/inner_product_computer.cc
@@ -51,16 +51,9 @@
   auto matrix = std::make_unique<CompressedRowSparseMatrix>(
       m_.num_cols(), m_.num_cols(), num_nonzeros);
   matrix->set_storage_type(storage_type);
-
   const CompressedRowBlockStructure* bs = m_.block_structure();
-  const std::vector<Block>& blocks = bs->cols;
-  matrix->mutable_row_blocks()->resize(blocks.size());
-  matrix->mutable_col_blocks()->resize(blocks.size());
-  for (int i = 0; i < blocks.size(); ++i) {
-    (*(matrix->mutable_row_blocks()))[i] = blocks[i].size;
-    (*(matrix->mutable_col_blocks()))[i] = blocks[i].size;
-  }
-
+  *matrix->mutable_row_blocks() = bs->cols;
+  *matrix->mutable_col_blocks() = bs->cols;
   return matrix;
 }
 
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 631862c..1825559 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -53,12 +53,6 @@
 #include "ceres/wall_time.h"
 
 namespace ceres::internal {
-
-using std::make_pair;
-using std::pair;
-using std::set;
-using std::vector;
-
 namespace {
 
 class BlockRandomAccessSparseMatrixAdapter final
@@ -174,12 +168,7 @@
     const CompressedRowBlockStructure* bs) {
   const int num_eliminate_blocks = options().elimination_groups[0];
   const int num_col_blocks = bs->cols.size();
-
-  vector<int> blocks(num_col_blocks - num_eliminate_blocks, 0);
-  for (int i = num_eliminate_blocks, j = 0; i < num_col_blocks; ++i, ++j) {
-    blocks[j] = bs->cols[i].size;
-  }
-
+  auto blocks = Tail(bs->cols, num_col_blocks - num_eliminate_blocks);
   set_lhs(std::make_unique<BlockRandomAccessDenseMatrix>(blocks));
   ResizeRhs(lhs()->num_rows());
 }
@@ -234,14 +223,11 @@
   const int num_col_blocks = bs->cols.size();
   const int num_row_blocks = bs->rows.size();
 
-  blocks_.resize(num_col_blocks - num_eliminate_blocks, 0);
-  for (int i = num_eliminate_blocks; i < num_col_blocks; ++i) {
-    blocks_[i - num_eliminate_blocks] = bs->cols[i].size;
-  }
+  blocks_ = Tail(bs->cols, num_col_blocks - num_eliminate_blocks);
 
-  set<pair<int, int>> block_pairs;
+  std::set<std::pair<int, int>> block_pairs;
   for (int i = 0; i < blocks_.size(); ++i) {
-    block_pairs.insert(make_pair(i, i));
+    block_pairs.emplace(i, i);
   }
 
   int r = 0;
@@ -250,7 +236,7 @@
     if (e_block_id >= num_eliminate_blocks) {
       break;
     }
-    vector<int> f_blocks;
+    std::vector<int> f_blocks;
 
     // Add to the chunk until the first block in the row is
     // different than the one in the first row for the chunk.
@@ -272,7 +258,7 @@
     f_blocks.erase(unique(f_blocks.begin(), f_blocks.end()), f_blocks.end());
     for (int i = 0; i < f_blocks.size(); ++i) {
       for (int j = i + 1; j < f_blocks.size(); ++j) {
-        block_pairs.insert(make_pair(f_blocks[i], f_blocks[j]));
+        block_pairs.emplace(f_blocks[i], f_blocks[j]);
       }
     }
   }
@@ -287,7 +273,7 @@
       for (const auto& cell : row.cells) {
         int r_block2_id = cell.block_id - num_eliminate_blocks;
         if (r_block1_id <= r_block2_id) {
-          block_pairs.insert(make_pair(r_block1_id, r_block2_id));
+          block_pairs.emplace(r_block1_id, r_block2_id);
         }
       }
     }
@@ -367,7 +353,7 @@
   // Extract block diagonal from the Schur complement to construct the
   // schur_jacobi preconditioner.
   for (int i = 0; i < blocks_.size(); ++i) {
-    const int block_size = blocks_[i];
+    const int block_size = blocks_[i].size;
 
     int sc_r, sc_c, sc_row_stride, sc_col_stride;
     CellInfo* sc_cell_info =
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h
index cac96b5..bb2cd88 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/internal/ceres/schur_complement_solver.h
@@ -183,8 +183,7 @@
   LinearSolver::Summary SolveReducedLinearSystemUsingConjugateGradients(
       const LinearSolver::PerSolveOptions& per_solve_options, double* solution);
 
-  // Size of the blocks in the Schur complement.
-  std::vector<int> blocks_;
+  std::vector<Block> blocks_;
   std::unique_ptr<SparseCholesky> sparse_cholesky_;
   std::unique_ptr<BlockRandomAccessDiagonalMatrix> preconditioner_;
   Vector cg_solution_;
diff --git a/internal/ceres/schur_eliminator_benchmark.cc b/internal/ceres/schur_eliminator_benchmark.cc
index efd33ec..271f92a 100644
--- a/internal/ceres/schur_eliminator_benchmark.cc
+++ b/internal/ceres/schur_eliminator_benchmark.cc
@@ -99,7 +99,8 @@
     b_.resize(matrix_->num_rows());
     b_.setRandom();
 
-    std::vector<int> blocks(1, kFBlockSize);
+    std::vector<Block> blocks;
+    blocks.emplace_back(kFBlockSize, 0);
     lhs_ = std::make_unique<BlockRandomAccessDenseMatrix>(blocks);
     diagonal_.resize(matrix_->num_cols());
     diagonal_.setOnes();
diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc
index c5ba328..dd30c23 100644
--- a/internal/ceres/schur_eliminator_test.cc
+++ b/internal/ceres/schur_eliminator_test.cc
@@ -57,8 +57,7 @@
 class SchurEliminatorTest : public ::testing::Test {
  protected:
   void SetUpFromId(int id) {
-    std::unique_ptr<LinearLeastSquaresProblem> problem =
-        CreateLinearLeastSquaresProblemFromId(id);
+    auto problem = CreateLinearLeastSquaresProblemFromId(id);
     CHECK(problem != nullptr);
     SetupHelper(problem.get());
   }
@@ -126,11 +125,7 @@
                                 const double relative_tolerance) {
     const CompressedRowBlockStructure* bs = A->block_structure();
     const int num_col_blocks = bs->cols.size();
-    std::vector<int> blocks(num_col_blocks - num_eliminate_blocks, 0);
-    for (int i = num_eliminate_blocks; i < num_col_blocks; ++i) {
-      blocks[i - num_eliminate_blocks] = bs->cols[i].size;
-    }
-
+    auto blocks = Tail(bs->cols, num_col_blocks - num_eliminate_blocks);
     BlockRandomAccessDenseMatrix lhs(blocks);
 
     const int num_cols = A->num_cols();
@@ -296,7 +291,8 @@
   Vector diagonal(matrix.num_cols());
   diagonal.setOnes();
 
-  std::vector<int> blocks(1, kFBlockSize);
+  std::vector<Block> blocks;
+  blocks.emplace_back(kFBlockSize, 0);
   BlockRandomAccessDenseMatrix actual_lhs(blocks);
   BlockRandomAccessDenseMatrix expected_lhs(blocks);
   Vector actual_rhs(kFBlockSize);
diff --git a/internal/ceres/schur_jacobi_preconditioner.cc b/internal/ceres/schur_jacobi_preconditioner.cc
index 3317927..9870924 100644
--- a/internal/ceres/schur_jacobi_preconditioner.cc
+++ b/internal/ceres/schur_jacobi_preconditioner.cc
@@ -51,9 +51,12 @@
                           << "SCHUR_JACOBI preconditioner.";
   CHECK(options_.context != nullptr);
 
-  std::vector<int> blocks(num_blocks);
+  std::vector<Block> blocks(num_blocks);
+  int position = 0;
   for (int i = 0; i < num_blocks; ++i) {
-    blocks[i] = bs.cols[i + options_.elimination_groups[0]].size;
+    blocks[i] =
+        Block(bs.cols[i + options_.elimination_groups[0]].size, position);
+    position += blocks[i].size;
   }
 
   m_ = std::make_unique<BlockRandomAccessDiagonalMatrix>(blocks);
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index 3f091b9..8e12874 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -61,9 +61,6 @@
 }
 }  // namespace
 
-using std::string;
-using std::vector;
-
 SuiteSparse::SuiteSparse() { cholmod_start(&cc_); }
 
 SuiteSparse::~SuiteSparse() { cholmod_finish(&cc_); }
@@ -163,7 +160,7 @@
 
 cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A,
                                              OrderingType ordering_type,
-                                             string* message) {
+                                             std::string* message) {
   cc_.nmethods = 1;
   cc_.method[0].ordering = OrderingTypeToCHOLMODEnum(ordering_type);
   cholmod_factor* factor = cholmod_analyze(A, &cc_);
@@ -183,7 +180,7 @@
 }
 
 cholmod_factor* SuiteSparse::AnalyzeCholeskyWithGivenOrdering(
-    cholmod_sparse* A, const vector<int>& ordering, string* message) {
+    cholmod_sparse* A, const std::vector<int>& ordering, std::string* message) {
   CHECK_EQ(ordering.size(), A->nrow);
 
   cc_.nmethods = 1;
@@ -207,9 +204,9 @@
 
 bool SuiteSparse::BlockOrdering(const cholmod_sparse* A,
                                 OrderingType ordering_type,
-                                const vector<int>& row_blocks,
-                                const vector<int>& col_blocks,
-                                vector<int>* ordering) {
+                                const std::vector<Block>& row_blocks,
+                                const std::vector<Block>& col_blocks,
+                                std::vector<int>* ordering) {
   if (ordering_type == OrderingType::NATURAL) {
     ordering->resize(A->nrow);
     for (int i = 0; i < A->nrow; ++i) {
@@ -223,8 +220,8 @@
 
   // Arrays storing the compressed column structure of the matrix
   // encoding the block sparsity of A.
-  vector<int> block_cols;
-  vector<int> block_rows;
+  std::vector<int> block_cols;
+  std::vector<int> block_rows;
 
   CompressedColumnScalarMatrixToBlockMatrix(reinterpret_cast<const int*>(A->i),
                                             reinterpret_cast<const int*>(A->p),
@@ -246,7 +243,7 @@
   block_matrix.sorted = 1;
   block_matrix.packed = 1;
 
-  vector<int> block_ordering(num_row_blocks);
+  std::vector<int> block_ordering(num_row_blocks);
   if (!Ordering(&block_matrix, ordering_type, block_ordering.data())) {
     return false;
   }
@@ -255,12 +252,13 @@
   return true;
 }
 
-cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(cholmod_sparse* A,
-                                                  OrderingType ordering_type,
-                                                  const vector<int>& row_blocks,
-                                                  const vector<int>& col_blocks,
-                                                  string* message) {
-  vector<int> ordering;
+cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(
+    cholmod_sparse* A,
+    OrderingType ordering_type,
+    const std::vector<Block>& row_blocks,
+    const std::vector<Block>& col_blocks,
+    std::string* message) {
+  std::vector<int> ordering;
   if (!BlockOrdering(A, ordering_type, row_blocks, col_blocks, &ordering)) {
     return nullptr;
   }
@@ -269,7 +267,7 @@
 
 LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A,
                                                   cholmod_factor* L,
-                                                  string* message) {
+                                                  std::string* message) {
   CHECK(A != nullptr);
   CHECK(L != nullptr);
 
@@ -328,7 +326,7 @@
 
 cholmod_dense* SuiteSparse::Solve(cholmod_factor* L,
                                   cholmod_dense* b,
-                                  string* message) {
+                                  std::string* message) {
   if (cc_.status != CHOLMOD_OK) {
     *message = "cholmod_solve failed. CHOLMOD status is not CHOLMOD_OK";
     return nullptr;
@@ -384,7 +382,7 @@
 }
 
 LinearSolverTerminationType SuiteSparseCholesky::Factorize(
-    CompressedRowSparseMatrix* lhs, string* message) {
+    CompressedRowSparseMatrix* lhs, std::string* message) {
   if (lhs == nullptr) {
     *message = "Failure: Input lhs is nullptr.";
     return LinearSolverTerminationType::FATAL_ERROR;
@@ -434,7 +432,7 @@
 
 LinearSolverTerminationType SuiteSparseCholesky::Solve(const double* rhs,
                                                        double* solution,
-                                                       string* message) {
+                                                       std::string* message) {
   // Error checking
   if (factor_ == nullptr) {
     *message = "Solve called without a call to Factorize first.";
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 14dd2b1..aeea7b2 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -44,6 +44,7 @@
 #include <vector>
 
 #include "SuiteSparseQR.hpp"
+#include "ceres/block_structure.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/linear_solver.h"
 #include "ceres/sparse_cholesky.h"
@@ -137,8 +138,8 @@
   // Block oriented version of AnalyzeCholesky.
   cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A,
                                        OrderingType ordering_type,
-                                       const std::vector<int>& row_blocks,
-                                       const std::vector<int>& col_blocks,
+                                       const std::vector<Block>& row_blocks,
+                                       const std::vector<Block>& col_blocks,
                                        std::string* message);
 
   // If A is symmetric, then compute the symbolic Cholesky
@@ -202,8 +203,8 @@
   // matrix to the numerical factorization routines.
   bool BlockOrdering(const cholmod_sparse* A,
                      OrderingType ordering_type,
-                     const std::vector<int>& row_blocks,
-                     const std::vector<int>& col_blocks,
+                     const std::vector<Block>& row_blocks,
+                     const std::vector<Block>& col_blocks,
                      std::vector<int>* ordering);
 
   // Nested dissection is only available if SuiteSparse is compiled
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index bdbe3c4..eb926ee 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -52,12 +52,6 @@
 
 namespace ceres::internal {
 
-using std::make_pair;
-using std::pair;
-using std::set;
-using std::swap;
-using std::vector;
-
 // TODO(sameeragarwal): Currently these are magic weights for the
 // preconditioner construction. Move these higher up into the Options
 // struct and provide some guidelines for choosing them.
@@ -81,10 +75,7 @@
   CHECK(options_.context != nullptr);
 
   // Vector of camera block sizes
-  block_size_.resize(num_blocks_);
-  for (int i = 0; i < num_blocks_; ++i) {
-    block_size_[i] = bs.cols[i + options_.elimination_groups[0]].size;
-  }
+  blocks_ = Tail(bs.cols, bs.cols.size() - options_.elimination_groups[0]);
 
   const time_t start_time = time(nullptr);
   switch (options_.type) {
@@ -124,13 +115,13 @@
 // preconditioner matrix.
 void VisibilityBasedPreconditioner::ComputeClusterJacobiSparsity(
     const CompressedRowBlockStructure& bs) {
-  vector<set<int>> visibility;
+  std::vector<std::set<int>> visibility;
   ComputeVisibility(bs, options_.elimination_groups[0], &visibility);
   CHECK_EQ(num_blocks_, visibility.size());
   ClusterCameras(visibility);
   cluster_pairs_.clear();
   for (int i = 0; i < num_clusters_; ++i) {
-    cluster_pairs_.insert(make_pair(i, i));
+    cluster_pairs_.insert(std::make_pair(i, i));
   }
 }
 
@@ -142,7 +133,7 @@
 // of edges in this forest are the cluster pairs.
 void VisibilityBasedPreconditioner::ComputeClusterTridiagonalSparsity(
     const CompressedRowBlockStructure& bs) {
-  vector<set<int>> visibility;
+  std::vector<std::set<int>> visibility;
   ComputeVisibility(bs, options_.elimination_groups[0], &visibility);
   CHECK_EQ(num_blocks_, visibility.size());
   ClusterCameras(visibility);
@@ -151,7 +142,7 @@
   // edges are the number of 3D points/e_blocks visible in both the
   // clusters at the ends of the edge. Return an approximate degree-2
   // maximum spanning forest of this graph.
-  vector<set<int>> cluster_visibility;
+  std::vector<std::set<int>> cluster_visibility;
   ComputeClusterVisibility(visibility, &cluster_visibility);
   auto cluster_graph = CreateClusterGraph(cluster_visibility);
   CHECK(cluster_graph != nullptr);
@@ -164,8 +155,7 @@
 void VisibilityBasedPreconditioner::InitStorage(
     const CompressedRowBlockStructure& bs) {
   ComputeBlockPairsInPreconditioner(bs);
-  m_ = std::make_unique<BlockRandomAccessSparseMatrix>(block_size_,
-                                                       block_pairs_);
+  m_ = std::make_unique<BlockRandomAccessSparseMatrix>(blocks_, block_pairs_);
 }
 
 // Call the canonical views algorithm and cluster the cameras based on
@@ -175,14 +165,14 @@
 // The cluster_membership_ vector is updated to indicate cluster
 // memberships for each camera block.
 void VisibilityBasedPreconditioner::ClusterCameras(
-    const vector<set<int>>& visibility) {
+    const std::vector<std::set<int>>& visibility) {
   auto schur_complement_graph = CreateSchurComplementGraph(visibility);
   CHECK(schur_complement_graph != nullptr);
 
   std::unordered_map<int, int> membership;
 
   if (options_.visibility_clustering_type == CANONICAL_VIEWS) {
-    vector<int> centers;
+    std::vector<int> centers;
     CanonicalViewsClusteringOptions clustering_options;
     clustering_options.size_penalty_weight = kCanonicalViewsSizePenaltyWeight;
     clustering_options.similarity_penalty_weight =
@@ -228,7 +218,7 @@
     const CompressedRowBlockStructure& bs) {
   block_pairs_.clear();
   for (int i = 0; i < num_blocks_; ++i) {
-    block_pairs_.insert(make_pair(i, i));
+    block_pairs_.insert(std::make_pair(i, i));
   }
 
   int r = 0;
@@ -256,7 +246,7 @@
       break;
     }
 
-    set<int> f_blocks;
+    std::set<int> f_blocks;
     for (; r < num_row_blocks; ++r) {
       const CompressedRow& row = bs.rows[r];
       if (row.cells.front().block_id != e_block_id) {
@@ -295,7 +285,7 @@
         const int block2 = cell.block_id - num_eliminate_blocks;
         if (block1 <= block2) {
           if (IsBlockPairInPreconditioner(block1, block2)) {
-            block_pairs_.insert(make_pair(block1, block2));
+            block_pairs_.insert(std::make_pair(block1, block2));
           }
         }
       }
@@ -392,7 +382,7 @@
     // dominance. See Lemma 1 in "Visibility Based Preconditioning
     // For Bundle Adjustment".
     MatrixRef m(cell_info->values, row_stride, col_stride);
-    m.block(r, c, block_size_[block1], block_size_[block2]) *= 0.5;
+    m.block(r, c, blocks_[block1].size, blocks_[block2].size) *= 0.5;
   }
 }
 
@@ -441,9 +431,9 @@
   int cluster1 = cluster_membership_[block1];
   int cluster2 = cluster_membership_[block2];
   if (cluster1 > cluster2) {
-    swap(cluster1, cluster2);
+    std::swap(cluster1, cluster2);
   }
-  return (cluster_pairs_.count(make_pair(cluster1, cluster2)) > 0);
+  return (cluster_pairs_.count(std::make_pair(cluster1, cluster2)) > 0);
 }
 
 bool VisibilityBasedPreconditioner::IsBlockPairOffDiagonal(
@@ -455,7 +445,7 @@
 // each vertex.
 void VisibilityBasedPreconditioner::ForestToClusterPairs(
     const WeightedGraph<int>& forest,
-    std::unordered_set<pair<int, int>, pair_hash>* cluster_pairs) const {
+    std::unordered_set<std::pair<int, int>, pair_hash>* cluster_pairs) const {
   CHECK(cluster_pairs != nullptr);
   cluster_pairs->clear();
   const std::unordered_set<int>& vertices = forest.vertices();
@@ -464,11 +454,11 @@
   // Add all the cluster pairs corresponding to the edges in the
   // forest.
   for (const int cluster1 : vertices) {
-    cluster_pairs->insert(make_pair(cluster1, cluster1));
+    cluster_pairs->insert(std::make_pair(cluster1, cluster1));
     const std::unordered_set<int>& neighbors = forest.Neighbors(cluster1);
     for (const int cluster2 : neighbors) {
       if (cluster1 < cluster2) {
-        cluster_pairs->insert(make_pair(cluster1, cluster2));
+        cluster_pairs->insert(std::make_pair(cluster1, cluster2));
       }
     }
   }
@@ -478,8 +468,8 @@
 // of all its cameras. In other words, the set of points visible to
 // any camera in the cluster.
 void VisibilityBasedPreconditioner::ComputeClusterVisibility(
-    const vector<set<int>>& visibility,
-    vector<set<int>>* cluster_visibility) const {
+    const std::vector<std::set<int>>& visibility,
+    std::vector<std::set<int>>* cluster_visibility) const {
   CHECK(cluster_visibility != nullptr);
   cluster_visibility->resize(0);
   cluster_visibility->resize(num_clusters_);
@@ -495,7 +485,7 @@
 // vertices.
 std::unique_ptr<WeightedGraph<int>>
 VisibilityBasedPreconditioner::CreateClusterGraph(
-    const vector<set<int>>& cluster_visibility) const {
+    const std::vector<std::set<int>>& cluster_visibility) const {
   auto cluster_graph = std::make_unique<WeightedGraph<int>>();
 
   for (int i = 0; i < num_clusters_; ++i) {
@@ -503,15 +493,15 @@
   }
 
   for (int i = 0; i < num_clusters_; ++i) {
-    const set<int>& cluster_i = cluster_visibility[i];
+    const std::set<int>& cluster_i = cluster_visibility[i];
     for (int j = i + 1; j < num_clusters_; ++j) {
-      vector<int> intersection;
-      const set<int>& cluster_j = cluster_visibility[j];
-      set_intersection(cluster_i.begin(),
-                       cluster_i.end(),
-                       cluster_j.begin(),
-                       cluster_j.end(),
-                       back_inserter(intersection));
+      std::vector<int> intersection;
+      const std::set<int>& cluster_j = cluster_visibility[j];
+      std::set_intersection(cluster_i.begin(),
+                            cluster_i.end(),
+                            cluster_j.begin(),
+                            cluster_j.end(),
+                            back_inserter(intersection));
 
       if (intersection.size() > 0) {
         // Clusters interact strongly when they share a large number
@@ -536,7 +526,7 @@
 // of integers so that the cluster ids are in [0, num_clusters_).
 void VisibilityBasedPreconditioner::FlattenMembershipMap(
     const std::unordered_map<int, int>& membership_map,
-    vector<int>* membership_vector) const {
+    std::vector<int>* membership_vector) const {
   CHECK(membership_vector != nullptr);
   membership_vector->resize(0);
   membership_vector->resize(num_blocks_, -1);
diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h
index f27eb9e..84773f8 100644
--- a/internal/ceres/visibility_based_preconditioner.h
+++ b/internal/ceres/visibility_based_preconditioner.h
@@ -55,6 +55,7 @@
 #include <utility>
 #include <vector>
 
+#include "ceres/block_structure.h"
 #include "ceres/graph.h"
 #include "ceres/linear_solver.h"
 #include "ceres/pair_hash.h"
@@ -176,7 +177,7 @@
   int num_clusters_;
 
   // Sizes of the blocks in the schur complement.
-  std::vector<int> block_size_;
+  std::vector<Block> blocks_;
 
   // Mapping from cameras to clusters.
   std::vector<int> cluster_membership_;