From c9d2ec8a9f6ceaf6a797716a88e41d2a7f799dcb Mon Sep 17 00:00:00 2001 From: Joydeep Biswas Date: Fri, 5 Aug 2022 13:53:42 -0500 Subject: [PATCH] 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 --- internal/ceres/block_sparse_matrix.cc | 37 +++++ internal/ceres/block_sparse_matrix.h | 2 + internal/ceres/block_sparse_matrix_test.cc | 141 +++++++++++++++++++ internal/ceres/block_structure.h | 76 ++++++++++ internal/ceres/triplet_sparse_matrix.cc | 32 +++++ internal/ceres/triplet_sparse_matrix.h | 4 + internal/ceres/triplet_sparse_matrix_test.cc | 39 +++++ 7 files changed, 331 insertions(+) diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index 9ba054f2c..22224a309 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc @@ -36,6 +36,7 @@ #include #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::ScaleColumns(const double* scale) { } } +void BlockSparseMatrix::ToCRSMatrix(CRSMatrix* crs_matrix) const { + CHECK(crs_matrix != nullptr); + crs_matrix->num_rows = num_rows_; + crs_matrix->num_cols = num_cols_; + vector& rows = crs_matrix->rows; + vector& cols = crs_matrix->cols; + vector& 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& 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 987a43562..da6b641f7 100644 --- a/internal/ceres/block_sparse_matrix.h +++ b/internal/ceres/block_sparse_matrix.h @@ -37,6 +37,7 @@ #include #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 @@ class CERES_NO_EXPORT BlockSparseMatrix final : public SparseMatrix { 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 d8cf98775..7b64a7c23 100644 --- a/internal/ceres/block_sparse_matrix_test.cc +++ b/internal/ceres/block_sparse_matrix_test.cc @@ -34,6 +34,7 @@ #include #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 +void CheckVectorEq(const std::vector& a, const std::vector& b) { + EXPECT_EQ(a.size(), b.size()); + for (int i = 0; i < a.size(); ++i) { + EXPECT_EQ(a[i], b[i]); + } +} + +std::unique_ptr 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 m = + std::make_unique(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 m = + std::make_unique(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, CreateDiagonalMatrix) { } } +TEST(BlockSparseMatrix, ToDenseMatrix) { + { + std::unique_ptr 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 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 m = CreateTestMatrixFromId(0); + CRSMatrix m_crs; + m->ToCRSMatrix(&m_crs); + std::vector rows_expected = {0, 2, 4, 7, 10}; + std::vector cols_expected = {0, 1, 0, 1, 2, 3, 4, 2, 3, 4}; + std::vector 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 m = CreateTestMatrixFromId(1); + CRSMatrix m_crs; + m->ToCRSMatrix(&m_crs); + std::vector rows_expected = {0, 4, 8, 9}; + std::vector cols_expected = {0, 1, 3, 4, 0, 1, 3, 4, 2}; + std::vector 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 52b7b2f0c..2ba6b327a 100644 --- a/internal/ceres/block_structure.h +++ b/internal/ceres/block_structure.h @@ -82,6 +82,82 @@ struct CERES_NO_EXPORT CompressedList { 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 cols; std::vector rows; diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc index c1c7ddcf7..86a444e26 100644 --- a/internal/ceres/triplet_sparse_matrix.cc +++ b/internal/ceres/triplet_sparse_matrix.cc @@ -33,6 +33,8 @@ #include #include +#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::ScaleColumns(const double* scale) { } } +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 @@ void TripletSparseMatrix::ToTextFile(FILE* file) const { } } +std::unique_ptr +TripletSparseMatrix::CreateFromTextFile(FILE* file) { + CHECK(file != nullptr); + int num_rows = 0; + int num_cols = 0; + std::vector rows; + std::vector cols; + std::vector 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( + num_rows, num_cols, rows, cols, values); +} + std::unique_ptr 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 2cf429086..a2ba1f198 100644 --- a/internal/ceres/triplet_sparse_matrix.h +++ b/internal/ceres/triplet_sparse_matrix.h @@ -34,6 +34,7 @@ #include #include +#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 @@ class CERES_NO_EXPORT TripletSparseMatrix final : public SparseMatrix { 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 @@ class CERES_NO_EXPORT TripletSparseMatrix final : public SparseMatrix { static std::unique_ptr CreateRandomMatrix( const TripletSparseMatrix::RandomMatrixOptions& options); + // Load a triplet sparse matrix from a text file. + static std::unique_ptr 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 e96c994c7..577d95990 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 @@ -343,4 +344,42 @@ TEST(TripletSparseMatrix, Resize) { } } +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