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_;