Updates to sparse block matrix structures to support new
sparse linear solvers.
* Add methods to convert TripletSparseMatrix and BlockSparseMatrix to
CRSMatrix structure.
* Added tests for conversion of TripletSparseMatrix and BlockSparseMatrix
to CRSMatrix structure.
* Added documentation on the BlockSparseMatrix structure.
Change-Id: I020cfa91c301567ceeb39ff2064183c5d88c9ed5
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index 9ba054f..22224a3 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -36,6 +36,7 @@
#include <vector>
#include "ceres/block_structure.h"
+#include "ceres/crs_matrix.h"
#include "ceres/internal/eigen.h"
#include "ceres/random.h"
#include "ceres/small_blas.h"
@@ -165,6 +166,42 @@
}
}
+void BlockSparseMatrix::ToCRSMatrix(CRSMatrix* crs_matrix) const {
+ CHECK(crs_matrix != nullptr);
+ crs_matrix->num_rows = num_rows_;
+ crs_matrix->num_cols = num_cols_;
+ vector<int>& rows = crs_matrix->rows;
+ vector<int>& cols = crs_matrix->cols;
+ vector<double>& values = crs_matrix->values;
+ rows.clear();
+ cols.clear();
+ values.clear();
+ rows.reserve(num_nonzeros_ + 1);
+ cols.reserve(num_nonzeros_);
+ values.reserve(num_nonzeros_);
+
+ for (const auto& row_block : block_structure_->rows) {
+ int row_block_size = row_block.block.size;
+ const vector<Cell>& cells = row_block.cells;
+ const int row_start = row_block.block.position;
+ for (int r = 0; r < row_block_size; ++r) {
+ rows.push_back(values.size());
+ for (const auto& cell : cells) {
+ int col_block_id = cell.block_id;
+ int col_block_size = block_structure_->cols[col_block_id].size;
+ int col_start = block_structure_->cols[col_block_id].position;
+ const MatrixRef m(
+ values_.get() + cell.position, row_block_size, col_block_size);
+ for (int c = 0; c < col_block_size; ++c) {
+ cols.push_back(c + col_start);
+ values.push_back(m(r, c));
+ }
+ }
+ }
+ }
+ rows.push_back(values.size());
+}
+
void BlockSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const {
CHECK(dense_matrix != nullptr);
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index 987a435..da6b641 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -37,6 +37,7 @@
#include <memory>
#include "ceres/block_structure.h"
+#include "ceres/crs_matrix.h"
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/export.h"
@@ -74,6 +75,7 @@
void LeftMultiply(const double* x, double* y) const final;
void SquaredColumnNorm(double* x) const final;
void ScaleColumns(const double* scale) final;
+ void ToCRSMatrix(CRSMatrix* matrix) const;
void ToDenseMatrix(Matrix* dense_matrix) const final;
void ToTextFile(FILE* file) const final;
diff --git a/internal/ceres/block_sparse_matrix_test.cc b/internal/ceres/block_sparse_matrix_test.cc
index d8cf987..7b64a7c 100644
--- a/internal/ceres/block_sparse_matrix_test.cc
+++ b/internal/ceres/block_sparse_matrix_test.cc
@@ -34,6 +34,7 @@
#include <string>
#include "ceres/casts.h"
+#include "ceres/crs_matrix.h"
#include "ceres/internal/eigen.h"
#include "ceres/linear_least_squares_problems.h"
#include "ceres/triplet_sparse_matrix.h"
@@ -43,6 +44,92 @@
namespace ceres {
namespace internal {
+namespace {
+template<typename T>
+void CheckVectorEq(const std::vector<T>& a, const std::vector<T>& b) {
+ EXPECT_EQ(a.size(), b.size());
+ for (int i = 0; i < a.size(); ++i) {
+ EXPECT_EQ(a[i], b[i]);
+ }
+}
+
+std::unique_ptr<BlockSparseMatrix> CreateTestMatrixFromId(int id) {
+ if (id == 0) {
+ // Create the following block sparse matrix:
+ // [ 1 2 0 0 0 0 ]
+ // [ 3 4 0 0 0 0 ]
+ // [ 0 0 5 6 7 0 ]
+ // [ 0 0 8 9 10 0 ]
+ CompressedRowBlockStructure* bs = new CompressedRowBlockStructure;
+ bs->cols = {
+ // Block size 2, position 0.
+ Block(2, 0),
+ // Block size 3, position 2.
+ Block(3, 2),
+ // Block size 1, position 5.
+ Block(1, 5),
+ };
+ bs->rows = {
+ CompressedRow(1),
+ CompressedRow(1)
+ };
+ bs->rows[0].block = Block(2, 0);
+ bs->rows[0].cells = { Cell(0, 0) };
+
+ bs->rows[1].block = Block(2, 2);
+ bs->rows[1].cells = { Cell(1, 4) };
+ std::unique_ptr<BlockSparseMatrix> m =
+ std::make_unique<BlockSparseMatrix>(bs);
+ EXPECT_NE(m, nullptr);
+ EXPECT_EQ(m->num_rows(), 4);
+ EXPECT_EQ(m->num_cols(), 6);
+ EXPECT_EQ(m->num_nonzeros(), 10);
+ double* values = m->mutable_values();
+ for (int i = 0; i < 10; ++i) {
+ values[i] = i + 1;
+ }
+ return m;
+ } else if (id == 1) {
+ // Create the following block sparse matrix:
+ // [ 1 2 0 5 6 0 ]
+ // [ 3 4 0 7 8 0 ]
+ // [ 0 0 9 0 0 0 ]
+ CompressedRowBlockStructure* bs = new CompressedRowBlockStructure;
+ bs->cols = {
+ // Block size 2, position 0.
+ Block(2, 0),
+ // Block size 1, position 2.
+ Block(1, 2),
+ // Block size 2, position 3.
+ Block(2, 3),
+ // Block size 1, position 5.
+ Block(1, 5),
+ };
+ bs->rows = {
+ CompressedRow(2),
+ CompressedRow(1)
+ };
+ bs->rows[0].block = Block(2, 0);
+ bs->rows[0].cells = { Cell(0, 0), Cell(2, 4) };
+
+ bs->rows[1].block = Block(1, 2);
+ bs->rows[1].cells = { Cell(1, 8) };
+ std::unique_ptr<BlockSparseMatrix> m =
+ std::make_unique<BlockSparseMatrix>(bs);
+ EXPECT_NE(m, nullptr);
+ EXPECT_EQ(m->num_rows(), 3);
+ EXPECT_EQ(m->num_cols(), 6);
+ EXPECT_EQ(m->num_nonzeros(), 9);
+ double* values = m->mutable_values();
+ for (int i = 0; i < 9; ++i) {
+ values[i] = i + 1;
+ }
+ return m;
+ }
+ return nullptr;
+}
+} // namespace
+
class BlockSparseMatrixTest : public ::testing::Test {
protected:
void SetUp() final {
@@ -214,5 +301,59 @@
}
}
+TEST(BlockSparseMatrix, ToDenseMatrix) {
+ {
+ std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(0);
+ Matrix m_dense;
+ m->ToDenseMatrix(&m_dense);
+ EXPECT_EQ(m_dense.rows(), 4);
+ EXPECT_EQ(m_dense.cols(), 6);
+ Matrix m_expected(4, 6);
+ m_expected << 1, 2, 0, 0, 0, 0,
+ 3, 4, 0, 0, 0, 0,
+ 0, 0, 5, 6, 7, 0,
+ 0, 0, 8, 9, 10, 0;
+ EXPECT_EQ(m_dense, m_expected);
+ }
+
+ {
+ std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(1);
+ Matrix m_dense;
+ m->ToDenseMatrix(&m_dense);
+ EXPECT_EQ(m_dense.rows(), 3);
+ EXPECT_EQ(m_dense.cols(), 6);
+ Matrix m_expected(3, 6);
+ m_expected << 1, 2, 0, 5, 6, 0,
+ 3, 4, 0, 7, 8, 0,
+ 0, 0, 9, 0, 0, 0;
+ EXPECT_EQ(m_dense, m_expected);
+ }
+}
+
+TEST(BlockSparseMatrix, ToCRSMatrix) {
+ {
+ std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(0);
+ CRSMatrix m_crs;
+ m->ToCRSMatrix(&m_crs);
+ std::vector<int> rows_expected = {0, 2, 4, 7, 10};
+ std::vector<int> cols_expected = {0, 1, 0, 1, 2, 3, 4, 2, 3, 4};
+ std::vector<double> values_expected = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+ CheckVectorEq(rows_expected, m_crs.rows);
+ CheckVectorEq(cols_expected, m_crs.cols);
+ CheckVectorEq(values_expected, m_crs.values);
+ }
+ {
+ std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(1);
+ CRSMatrix m_crs;
+ m->ToCRSMatrix(&m_crs);
+ std::vector<int> rows_expected = {0, 4, 8, 9};
+ std::vector<int> cols_expected = {0, 1, 3, 4, 0, 1, 3, 4, 2};
+ std::vector<double> values_expected = {1, 2, 5, 6, 3, 4, 7, 8, 9};
+ CheckVectorEq(rows_expected, m_crs.rows);
+ CheckVectorEq(cols_expected, m_crs.cols);
+ CheckVectorEq(values_expected, m_crs.values);
+ }
+}
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/block_structure.h b/internal/ceres/block_structure.h
index 52b7b2f..2ba6b32 100644
--- a/internal/ceres/block_structure.h
+++ b/internal/ceres/block_structure.h
@@ -82,6 +82,82 @@
using CompressedRow = CompressedList;
using CompressedColumn = CompressedList;
+// CompressedRowBlockStructure specifies the storage structure of a row block
+// sparse matrix.
+//
+// Consider the following matrix A:
+// A = [A_11 A_12 ...
+// A_21 A_22 ...
+// ...
+// A_m1 A_m2 ... ]
+//
+// A row block sparse matrix is a matrix where the following properties hold:
+// 1. The number of rows in every block A_ij and A_ik are the same.
+// 2. The number of columns in every block A_ij and A_kj are the same.
+// 3. The number of rows in A_ij and A_kj may be different (i != k).
+// 4. The number of columns in A_ij and A_ik may be different (j != k).
+// 5. Any block A_ij may be all 0s, in which case the block is not stored.
+//
+// The structure of the matrix is stored as follows:
+//
+// The `rows' array contains the following information for each row block:
+// - rows[i].block.size: The number of rows in each block A_ij in the row block.
+// - rows[i].block.position: The starting row in the full matrix A of the
+// row block i.
+// - rows[i].cells[j].block_id: The index into the `cols' array corresponding to
+// the non-zero blocks A_ij.
+// - rows[i].cells[j].position: The index in the `values' array for the contents
+// of block A_ij.
+//
+// The `cols' array contains the following information for block:
+// - cols[.].size: The number of columns spanned by the block.
+// - cols[.].position: The starting column in the full matrix A of the block.
+//
+//
+// Example of a row block sparse matrix:
+// block_id: | 0 |1|2 |3 |
+// rows[0]: [ 1 2 0 3 4 0 ]
+// [ 5 6 0 7 8 0 ]
+// rows[1]: [ 0 0 9 0 0 0 ]
+//
+// This matrix is stored as follows:
+//
+// There are four column blocks:
+// cols[0].size = 2
+// cols[0].position = 0
+// cols[1].size = 1
+// cols[1].position = 2
+// cols[2].size = 2
+// cols[2].position = 3
+// cols[3].size = 1
+// cols[3].position = 5
+
+// The first row block spans two rows, starting at row 0:
+// rows[0].block.size = 2 // This row block spans two rows.
+// rows[0].block.position = 0 // It starts at row 0.
+// rows[0] has two cells, at column blocks 0 and 2:
+// rows[0].cells[0].block_id = 0 // This cell is in column block 0.
+// rows[0].cells[0].position = 0 // See below for an explanation of this.
+// rows[0].cells[1].block_id = 2 // This cell is in column block 2.
+// rows[0].cells[1].position = 4 // See below for an explanation of this.
+//
+// The second row block spans two rows, starting at row 2:
+// rows[1].block.size = 1 // This row block spans one row.
+// rows[1].block.position = 2 // It starts at row 2.
+// rows[1] has one cell at column block 1:
+// rows[1].cells[0].block_id = 1 // This cell is in column block 1.
+// rows[1].cells[0].position = 8 // See below for an explanation of this.
+//
+// The values in each blocks are stored contiguously in row major order.
+// However, there is no unique way to order the blocks -- it is usually
+// optimized to promote cache coherent access, e.g. ordering it so that
+// Jacobian blocks of parameters of the same type are stored nearby.
+// This is one possible way to store the values of the blocks in a values array:
+// values = { 1, 2, 5, 6, 3, 4, 7, 8, 9 }
+// | | | | // The three blocks.
+// ^ rows[0].cells[0].position = 0
+// ^ rows[0].cells[1].position = 4
+// ^ rows[1].cells[0].position = 8
struct CERES_NO_EXPORT CompressedRowBlockStructure {
std::vector<Block> cols;
std::vector<CompressedRow> rows;
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc
index c1c7ddc..86a444e 100644
--- a/internal/ceres/triplet_sparse_matrix.cc
+++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -33,6 +33,8 @@
#include <algorithm>
#include <memory>
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/crs_matrix.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/export.h"
#include "ceres/random.h"
@@ -194,6 +196,11 @@
}
}
+void TripletSparseMatrix::ToCRSMatrix(CRSMatrix* crs_matrix) const {
+ CompressedRowSparseMatrix::FromTripletSparseMatrix(*this)->ToCRSMatrix(
+ crs_matrix);
+}
+
void TripletSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const {
dense_matrix->resize(num_rows_, num_cols_);
dense_matrix->setZero();
@@ -275,6 +282,31 @@
}
}
+std::unique_ptr<TripletSparseMatrix>
+TripletSparseMatrix::CreateFromTextFile(FILE* file) {
+ CHECK(file != nullptr);
+ int num_rows = 0;
+ int num_cols = 0;
+ std::vector<int> rows;
+ std::vector<int> cols;
+ std::vector<double> values;
+ while (true) {
+ int row, col;
+ double value;
+ if (fscanf(file, "%d %d %lf", &row, &col, &value) != 3) {
+ break;
+ }
+ rows.push_back(row);
+ cols.push_back(col);
+ values.push_back(value);
+ num_rows = std::max(num_rows, row + 1);
+ num_cols = std::max(num_cols, col + 1);
+ }
+ VLOG(1) << "Read " << rows.size() << " nonzeros from file.";
+ return std::make_unique<TripletSparseMatrix>(
+ num_rows, num_cols, rows, cols, values);
+}
+
std::unique_ptr<TripletSparseMatrix> TripletSparseMatrix::CreateRandomMatrix(
const TripletSparseMatrix::RandomMatrixOptions& options) {
CHECK_GT(options.num_rows, 0);
diff --git a/internal/ceres/triplet_sparse_matrix.h b/internal/ceres/triplet_sparse_matrix.h
index 2cf4290..a2ba1f1 100644
--- a/internal/ceres/triplet_sparse_matrix.h
+++ b/internal/ceres/triplet_sparse_matrix.h
@@ -34,6 +34,7 @@
#include <memory>
#include <vector>
+#include "ceres/crs_matrix.h"
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/export.h"
@@ -68,6 +69,7 @@
void LeftMultiply(const double* x, double* y) const final;
void SquaredColumnNorm(double* x) const final;
void ScaleColumns(const double* scale) final;
+ void ToCRSMatrix(CRSMatrix* matrix) const;
void ToDenseMatrix(Matrix* dense_matrix) const final;
void ToTextFile(FILE* file) const final;
// clang-format off
@@ -135,6 +137,8 @@
static std::unique_ptr<TripletSparseMatrix> CreateRandomMatrix(
const TripletSparseMatrix::RandomMatrixOptions& options);
+ // Load a triplet sparse matrix from a text file.
+ static std::unique_ptr<TripletSparseMatrix> CreateFromTextFile(FILE* file);
private:
void AllocateMemory();
void CopyData(const TripletSparseMatrix& orig);
diff --git a/internal/ceres/triplet_sparse_matrix_test.cc b/internal/ceres/triplet_sparse_matrix_test.cc
index e96c994..577d959 100644
--- a/internal/ceres/triplet_sparse_matrix_test.cc
+++ b/internal/ceres/triplet_sparse_matrix_test.cc
@@ -28,6 +28,7 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
+#include "ceres/crs_matrix.h"
#include "ceres/triplet_sparse_matrix.h"
#include <memory>
@@ -343,4 +344,42 @@
}
}
+TEST(TripletSparseMatrix, ToCRSMatrix) {
+ // Test matrix:
+ // [1, 2, 0, 5, 6, 0,
+ // 3, 4, 0, 7, 8, 0,
+ // 0, 0, 9, 0, 0, 0]
+ TripletSparseMatrix m(3,
+ 6,
+ {0, 0, 0, 0, 1, 1, 1, 1, 2},
+ {0, 1, 3, 4, 0, 1, 3, 4, 2},
+ {1, 2, 3, 4, 5, 6, 7, 8, 9});
+ CRSMatrix m_crs;
+ m.ToCRSMatrix(&m_crs);
+ EXPECT_EQ(m_crs.num_rows, 3);
+ EXPECT_EQ(m_crs.num_cols, 6);
+
+ EXPECT_EQ(m_crs.rows.size(), 4);
+ EXPECT_EQ(m_crs.rows[0], 0);
+ EXPECT_EQ(m_crs.rows[1], 4);
+ EXPECT_EQ(m_crs.rows[2], 8);
+ EXPECT_EQ(m_crs.rows[3], 9);
+
+ EXPECT_EQ(m_crs.cols.size(), 9);
+ EXPECT_EQ(m_crs.cols[0], 0);
+ EXPECT_EQ(m_crs.cols[1], 1);
+ EXPECT_EQ(m_crs.cols[2], 3);
+ EXPECT_EQ(m_crs.cols[3], 4);
+ EXPECT_EQ(m_crs.cols[4], 0);
+ EXPECT_EQ(m_crs.cols[5], 1);
+ EXPECT_EQ(m_crs.cols[6], 3);
+ EXPECT_EQ(m_crs.cols[7], 4);
+ EXPECT_EQ(m_crs.cols[8], 2);
+
+ EXPECT_EQ(m_crs.values.size(), 9);
+ for (int i = 0; i < 9; ++i) {
+ EXPECT_EQ(m_crs.values[i], i + 1);
+ }
+}
+
} // namespace ceres::internal
