From 1ece5a95fbb4d8f089d5da30eff4e9d831db42fe Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Thu, 22 Jun 2017 00:21:09 -0700 Subject: [PATCH] Delete obsolete code Remove outer product computation code from CompressedRowSparseMatrix. In the process also remove the crsb_cols and crsb_rows vectors from the matrix, which were added to carry the block sparsity of the matrix so that the outer product could be computed fast. InnerProductComputer and its reliance on BlockSparseMatrix has rendered all of this code moot. Change-Id: If3ee0dc8ad4ff79594fd1eebc15a647c4495d726 --- .../ceres/compressed_row_jacobian_writer.cc | 11 - .../ceres/compressed_row_sparse_matrix.cc | 422 ------------------ internal/ceres/compressed_row_sparse_matrix.h | 52 --- .../compressed_row_sparse_matrix_test.cc | 124 ----- 4 files changed, 609 deletions(-) diff --git a/internal/ceres/compressed_row_jacobian_writer.cc b/internal/ceres/compressed_row_jacobian_writer.cc index d21aa7d65..0444d45fe 100644 --- a/internal/ceres/compressed_row_jacobian_writer.cc +++ b/internal/ceres/compressed_row_jacobian_writer.cc @@ -123,12 +123,6 @@ SparseMatrix* CompressedRowJacobianWriter::CreateJacobian() const { int* rows = jacobian->mutable_rows(); int* cols = jacobian->mutable_cols(); - // Initialize crsb rows and cols. - std::vector& crsb_rows = *jacobian->mutable_crsb_rows(); - std::vector& crsb_cols = *jacobian->mutable_crsb_cols(); - crsb_rows.resize(residual_blocks.size() + 1); - crsb_rows[0] = 0; - int row_pos = 0; rows[0] = 0; for (int i = 0; i < residual_blocks.size(); ++i) { @@ -165,11 +159,6 @@ SparseMatrix* CompressedRowJacobianWriter::CreateJacobian() const { << "Parameter Blocks: " << parameter_block_description; } - // Populate crsb rows and cols. - crsb_rows[i + 1] = crsb_rows[i] + parameter_indices.size(); - std::copy(parameter_indices.begin(), parameter_indices.end(), - std::back_inserter(crsb_cols)); - // Update the row indices. const int num_residuals = residual_block->NumResiduals(); for (int j = 0; j < num_residuals; ++j) { diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc index 3e75f1922..ee59eae43 100644 --- a/internal/ceres/compressed_row_sparse_matrix.cc +++ b/internal/ceres/compressed_row_sparse_matrix.cc @@ -311,10 +311,6 @@ void CompressedRowSparseMatrix::DeleteRows(int delta_rows) { return; } - // Sanity check for compressed row sparse block information - CHECK_EQ(crsb_rows_.size(), row_blocks_.size() + 1); - CHECK_EQ(crsb_rows_.back(), crsb_cols_.size()); - // Walk the list of row blocks until we reach the new number of rows // and the drop the rest of the row blocks. int num_row_blocks = 0; @@ -325,11 +321,6 @@ void CompressedRowSparseMatrix::DeleteRows(int delta_rows) { } row_blocks_.resize(num_row_blocks); - - // Update compressed row sparse block (crsb) information. - CHECK_EQ(num_rows, num_rows_); - crsb_rows_.resize(num_row_blocks + 1); - crsb_cols_.resize(crsb_rows_[num_row_blocks]); } void CompressedRowSparseMatrix::AppendRows(const CompressedRowSparseMatrix& m) { @@ -377,32 +368,8 @@ void CompressedRowSparseMatrix::AppendRows(const CompressedRowSparseMatrix& m) { return; } - // Sanity check for compressed row sparse block information - CHECK_EQ(crsb_rows_.size(), row_blocks_.size() + 1); - CHECK_EQ(crsb_rows_.back(), crsb_cols_.size()); - row_blocks_.insert( row_blocks_.end(), m.row_blocks().begin(), m.row_blocks().end()); - - // The rest of the code updates the compressed row sparse block - // (crsb) information. - const int num_crsb_nonzeros = crsb_cols_.size(); - const int m_num_crsb_nonzeros = m.crsb_cols_.size(); - crsb_cols_.resize(num_crsb_nonzeros + m_num_crsb_nonzeros); - std::copy(&m.crsb_cols()[0], - &m.crsb_cols()[0] + m_num_crsb_nonzeros, - &crsb_cols_[num_crsb_nonzeros]); - - const int num_crsb_rows = crsb_rows_.size() - 1; - const int m_num_crsb_rows = m.crsb_rows_.size() - 1; - crsb_rows_.resize(num_crsb_rows + m_num_crsb_rows + 1); - std::fill(crsb_rows_.begin() + num_crsb_rows, - crsb_rows_.begin() + num_crsb_rows + m_num_crsb_rows + 1, - crsb_rows_[num_crsb_rows]); - - for (int r = 0; r < m_num_crsb_rows + 1; ++r) { - crsb_rows_[num_crsb_rows + r] += m.crsb_rows()[r]; - } } void CompressedRowSparseMatrix::ToTextFile(FILE* file) const { @@ -469,15 +436,6 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateBlockDiagonalMatrix( *matrix->mutable_row_blocks() = blocks; *matrix->mutable_col_blocks() = blocks; - // Fill compressed row sparse block (crsb) information. - vector& crsb_rows = *matrix->mutable_crsb_rows(); - vector& crsb_cols = *matrix->mutable_crsb_cols(); - for (int i = 0; i < blocks.size(); ++i) { - crsb_rows.push_back(i); - crsb_cols.push_back(i); - } - crsb_rows.push_back(blocks.size()); - CHECK_EQ(idx_cursor, num_nonzeros); CHECK_EQ(col_cursor, num_rows); return matrix; @@ -517,381 +475,11 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::Transpose() const { return transpose; } - // Sanity check for compressed row sparse block information - CHECK_EQ(crsb_rows_.size(), row_blocks_.size() + 1); - CHECK_EQ(crsb_rows_.back(), crsb_cols_.size()); - *(transpose->mutable_row_blocks()) = col_blocks_; *(transpose->mutable_col_blocks()) = row_blocks_; - - // The rest of the code updates the compressed row sparse block - // (crsb) information. - vector& transpose_crsb_rows = *transpose->mutable_crsb_rows(); - vector& transpose_crsb_cols = *transpose->mutable_crsb_cols(); - - transpose_crsb_rows.resize(col_blocks_.size() + 1); - transpose_crsb_cols.resize(crsb_cols_.size()); - TransposeForCompressedRowSparseStructure(row_blocks().size(), - col_blocks().size(), - crsb_cols().size(), - crsb_rows().data(), - crsb_cols().data(), - NULL, - transpose_crsb_rows.data(), - transpose_crsb_cols.data(), - NULL); - return transpose; } -namespace { -// A ProductTerm is a term in the block outer product of a matrix with -// itself. -struct ProductTerm { - ProductTerm(const int row, const int col, const int index) - : row(row), col(col), index(index) {} - - bool operator<(const ProductTerm& right) const { - if (row == right.row) { - if (col == right.col) { - return index < right.index; - } - return col < right.col; - } - return row < right.row; - } - - int row; - int col; - int index; -}; - -// Create outer product matrix based on the block product information. -// The input block product is already sorted. This function does not -// set the sparse rows/cols information. Instead, it only collects the -// nonzeros for each compressed row and puts in row_nnz. The caller of -// this function will traverse the block product in a second round to -// generate the sparse rows/cols information. This function also -// computes the block offset information for the outer product matrix, -// which is used in outer product computation. -CompressedRowSparseMatrix* CreateOuterProductMatrix( - const int num_cols, - const CompressedRowSparseMatrix::StorageType storage_type, - const vector& blocks, - const vector& product, - vector* row_nnz) { - // Count the number of unique product term, which in turn is the - // number of non-zeros in the outer product. Also count the number - // of non-zeros in each row. - row_nnz->resize(blocks.size()); - std::fill(row_nnz->begin(), row_nnz->end(), 0); - (*row_nnz)[product[0].row] = blocks[product[0].col]; - int num_nonzeros = blocks[product[0].row] * blocks[product[0].col]; - for (int i = 1; i < product.size(); ++i) { - // Each (row, col) block counts only once. - // This check depends on product sorted on (row, col). - if (product[i].row != product[i - 1].row || - product[i].col != product[i - 1].col) { - (*row_nnz)[product[i].row] += blocks[product[i].col]; - num_nonzeros += blocks[product[i].row] * blocks[product[i].col]; - } - } - - CompressedRowSparseMatrix* matrix = - new CompressedRowSparseMatrix(num_cols, num_cols, num_nonzeros); - matrix->set_storage_type(storage_type); - - *(matrix->mutable_row_blocks()) = blocks; - *(matrix->mutable_col_blocks()) = blocks; - - // Compute block offsets for outer product matrix, which is used in - // ComputeOuterProduct. - vector* block_offsets = matrix->mutable_block_offsets(); - block_offsets->resize(blocks.size() + 1); - (*block_offsets)[0] = 0; - for (int i = 0; i < blocks.size(); ++i) { - (*block_offsets)[i + 1] = (*block_offsets)[i] + blocks[i]; - } - - return matrix; -} - -CompressedRowSparseMatrix* CompressAndFillProgram( - const int num_cols, - const CompressedRowSparseMatrix::StorageType storage_type, - const vector& blocks, - const vector& product, - vector* program) { - CHECK_GT(product.size(), 0); - - vector row_nnz; - CompressedRowSparseMatrix* matrix = CreateOuterProductMatrix( - num_cols, storage_type, blocks, product, &row_nnz); - - const vector& block_offsets = matrix->block_offsets(); - - int* crsm_rows = matrix->mutable_rows(); - std::fill(crsm_rows, crsm_rows + num_cols + 1, 0); - int* crsm_cols = matrix->mutable_cols(); - std::fill(crsm_cols, crsm_cols + matrix->num_nonzeros(), 0); - - CHECK_NOTNULL(program)->clear(); - program->resize(product.size()); - - // Non zero elements are not stored consecutively across rows in a block. - // We seperate nonzero into three categories: - // nonzeros in all previous row blocks counted in nnz - // nonzeros in current row counted in row_nnz - // nonzeros in previous col blocks of current row counted in col_nnz - // - // Give an element (j, k) within a block such that j and k - // represent the relative position to the starting row and starting col of - // the block, the row and col for the element is - // block_offsets[current.row] + j - // block_offsets[current.col] + k - // The total number of nonzero to the element is - // nnz + row_nnz[current.row] * j + col_nnz + k - // - // program keeps col_nnz for block product, which is used later for - // outer product computation. - // - // There is no special handling for diagonal blocks as we generate - // BLOCK triangular matrix (diagonal block is full block) instead of - // standard triangular matrix. - int nnz = 0; - int col_nnz = 0; - - // Process first product term. - for (int j = 0; j < blocks[product[0].row]; ++j) { - crsm_rows[block_offsets[product[0].row] + j + 1] = row_nnz[product[0].row]; - for (int k = 0; k < blocks[product[0].col]; ++k) { - crsm_cols[row_nnz[product[0].row] * j + k] = - block_offsets[product[0].col] + k; - } - } - - (*program)[product[0].index] = 0; - - // Process rest product terms. - for (int i = 1; i < product.size(); ++i) { - const ProductTerm& previous = product[i - 1]; - const ProductTerm& current = product[i]; - - // Sparsity structure is updated only if the term is not a repeat. - if (previous.row != current.row || previous.col != current.col) { - col_nnz += blocks[previous.col]; - if (previous.row != current.row) { - nnz += col_nnz * blocks[previous.row]; - col_nnz = 0; - - for (int j = 0; j < blocks[current.row]; ++j) { - crsm_rows[block_offsets[current.row] + j + 1] = row_nnz[current.row]; - } - } - - for (int j = 0; j < blocks[current.row]; ++j) { - for (int k = 0; k < blocks[current.col]; ++k) { - crsm_cols[nnz + row_nnz[current.row] * j + col_nnz + k] = - block_offsets[current.col] + k; - } - } - } - - (*program)[current.index] = col_nnz; - } - - for (int i = 1; i < num_cols + 1; ++i) { - crsm_rows[i] += crsm_rows[i - 1]; - } - - return matrix; -} - -// input is a matrix of dimesion -// output is a matrix of dimension -// -// Implement block multiplication O = I1' * I2. -// I1 is block(0, col_block1_begin, row_block_size, col_block1_size) of input -// I2 is block(0, col_block2_begin, row_block_size, col_block2_size) of input -// O is block(0, 0, col_block1_size, col_block2_size) of output -void ComputeBlockMultiplication(const int row_block_size, - const int col_block1_size, - const int col_block2_size, - const int col_block1_begin, - const int col_block2_begin, - const int input_cols, - const double* input, - const int output_cols, - double* output) { - for (int r = 0; r < row_block_size; ++r) { - for (int idx1 = 0; idx1 < col_block1_size; ++idx1) { - for (int idx2 = 0; idx2 < col_block2_size; ++idx2) { - output[output_cols * idx1 + idx2] += - input[input_cols * r + col_block1_begin + idx1] * - input[input_cols * r + col_block2_begin + idx2]; - } - } - } -} -} // namespace - -CompressedRowSparseMatrix* -CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( - const CompressedRowSparseMatrix& m, - const CompressedRowSparseMatrix::StorageType storage_type, - vector* program) { - CHECK(storage_type == LOWER_TRIANGULAR || storage_type == UPPER_TRIANGULAR); - CHECK_NOTNULL(program)->clear(); - CHECK_GT(m.num_nonzeros(), 0) - << "Congratulations, you found a bug in Ceres. Please report it."; - - vector product; - const vector& col_blocks = m.col_blocks(); - const vector& crsb_rows = m.crsb_rows(); - const vector& crsb_cols = m.crsb_cols(); - - // Give input matrix m in Compressed Row Sparse Block format - // (row_block, col_block) - // represent each block multiplication - // (row_block, col_block1)' X (row_block, col_block2) - // by its product term index and sort the product terms - // (col_block1, col_block2, index) - // - // Due to the compression on rows, col_block is accessed through idx to - // crsb_cols. So col_block is accessed as crsb_cols[idx] in the code. - for (int row_block = 1; row_block < crsb_rows.size(); ++row_block) { - for (int idx1 = crsb_rows[row_block - 1]; idx1 < crsb_rows[row_block]; - ++idx1) { - if (storage_type == LOWER_TRIANGULAR) { - for (int idx2 = crsb_rows[row_block - 1]; idx2 <= idx1; ++idx2) { - product.push_back( - ProductTerm(crsb_cols[idx1], crsb_cols[idx2], product.size())); - } - } else { // Upper triangular matrix. - for (int idx2 = idx1; idx2 < crsb_rows[row_block]; ++idx2) { - product.push_back( - ProductTerm(crsb_cols[idx1], crsb_cols[idx2], product.size())); - } - } - } - } - - sort(product.begin(), product.end()); - return CompressAndFillProgram( - m.num_cols(), storage_type, col_blocks, product, program); -} - -// Give input matrix m in Compressed Row Sparse Block format -// (row_block, col_block) -// compute outer product m' * m as sum of block multiplications -// (row_block, col_block1)' X (row_block, col_block2) -// -// Given row_block of the input matrix m, we use m_row_begin to represent -// the starting row of the row block and m_row_nnz to represent number of -// nonzero in each row of the row block, then the rows belonging to -// the row block can be represented as a dense matrix starting at -// m.values() + m.rows()[m_row_begin] -// with dimension -// -// -// Then each input matrix block (row_block, col_block) can be represented as -// a block of above dense matrix starting at position -// (0, m_col_nnz) -// with size -// -// where m_col_nnz is the number of nonzero before col_block in each row. -// -// The outer product block is represented similarly with m_row_begin, -// m_row_nnz, m_col_nnz, etc. replaced by row_begin, row_nnz, col_nnz, -// etc. The difference is, m_row_begin and m_col_nnz is counted -// during the traverse of block multiplication, while row_begin and -// col_nnz are got from pre-computed block_offsets and program. -// -// Due to the compression on rows, col_block is accessed through -// idx to crsb_col vector. So col_block is accessed as crsb_col[idx] -// in the code. -// -// Note this function produces a triangular matrix in block unit (i.e. -// diagonal block is a normal block) instead of standard triangular matrix. -// So there is no special handling for diagonal blocks. -void CompressedRowSparseMatrix::ComputeOuterProduct( - const CompressedRowSparseMatrix& m, - const vector& program, - CompressedRowSparseMatrix* result) { - CHECK(result->storage_type() == LOWER_TRIANGULAR || - result->storage_type() == UPPER_TRIANGULAR); - result->SetZero(); - double* values = result->mutable_values(); - const int* rows = result->rows(); - const vector& block_offsets = result->block_offsets(); - - int cursor = 0; - const double* m_values = m.values(); - const int* m_rows = m.rows(); - const vector& row_blocks = m.row_blocks(); - const vector& col_blocks = m.col_blocks(); - const vector& crsb_rows = m.crsb_rows(); - const vector& crsb_cols = m.crsb_cols(); - const StorageType storage_type = result->storage_type(); -#define COL_BLOCK1 (crsb_cols[idx1]) -#define COL_BLOCK2 (crsb_cols[idx2]) - - // Iterate row blocks. - for (int row_block = 0, m_row_begin = 0; row_block < row_blocks.size(); - m_row_begin += row_blocks[row_block++]) { - // Non zeros are not stored consecutively across rows in a block. - // The gaps between rows is the number of nonzeros of the - // input matrix compressed row. - const int m_row_nnz = m_rows[m_row_begin + 1] - m_rows[m_row_begin]; - - // Iterate (col_block1 x col_block2). - for (int idx1 = crsb_rows[row_block], m_col_nnz1 = 0; - idx1 < crsb_rows[row_block + 1]; - m_col_nnz1 += col_blocks[COL_BLOCK1], ++idx1) { - // Non zeros are not stored consecutively across rows in a - // block. The gaps between rows is the number of nonzeros of the - // outer product matrix compressed row. - const int row_begin = block_offsets[COL_BLOCK1]; - const int row_nnz = rows[row_begin + 1] - rows[row_begin]; - if (storage_type == LOWER_TRIANGULAR) { - for (int idx2 = crsb_rows[row_block], m_col_nnz2 = 0; idx2 <= idx1; - m_col_nnz2 += col_blocks[COL_BLOCK2], ++idx2, ++cursor) { - int col_nnz = program[cursor]; - ComputeBlockMultiplication(row_blocks[row_block], - col_blocks[COL_BLOCK1], - col_blocks[COL_BLOCK2], - m_col_nnz1, - m_col_nnz2, - m_row_nnz, - m_values + m_rows[m_row_begin], - row_nnz, - values + rows[row_begin] + col_nnz); - } - } else { - for (int idx2 = idx1, m_col_nnz2 = m_col_nnz1; - idx2 < crsb_rows[row_block + 1]; - m_col_nnz2 += col_blocks[COL_BLOCK2], ++idx2, ++cursor) { - int col_nnz = program[cursor]; - ComputeBlockMultiplication(row_blocks[row_block], - col_blocks[COL_BLOCK1], - col_blocks[COL_BLOCK2], - m_col_nnz1, - m_col_nnz2, - m_row_nnz, - m_values + m_rows[m_row_begin], - row_nnz, - values + rows[row_begin] + col_nnz); - } - } - } - } - -#undef COL_BLOCK1 -#undef COL_BLOCK2 - - CHECK_EQ(cursor, program.size()); -} - CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix( const CompressedRowSparseMatrix::RandomMatrixOptions& options) { CHECK_GT(options.num_row_blocks, 0); @@ -924,8 +512,6 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix( col_blocks.push_back(options.min_col_block_size + delta_block_size); } - vector crsb_rows; - vector crsb_cols; vector tsm_rows; vector tsm_cols; vector tsm_values; @@ -939,8 +525,6 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix( // not what the user wants, so do the matrix generation till we have // at least one non-zero entry. while (tsm_values.empty()) { - crsb_rows.clear(); - crsb_cols.clear(); tsm_rows.clear(); tsm_cols.clear(); tsm_values.clear(); @@ -948,7 +532,6 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix( int row_block_begin = 0; for (int r = 0; r < options.num_row_blocks; ++r) { int col_block_begin = 0; - crsb_rows.push_back(crsb_cols.size()); for (int c = 0; c < options.num_col_blocks; ++c) { // Randomly determine if this block is present or not. if (RandDouble() <= options.block_density) { @@ -959,14 +542,11 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix( &tsm_rows, &tsm_cols, &tsm_values); - // Add the block to the block sparse structure. - crsb_cols.push_back(c); } col_block_begin += col_blocks[c]; } row_block_begin += row_blocks[r]; } - crsb_rows.push_back(crsb_cols.size()); } const int num_rows = std::accumulate(row_blocks.begin(), row_blocks.end(), 0); @@ -979,8 +559,6 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix( kDoNotTranspose); (*matrix->mutable_row_blocks()) = row_blocks; (*matrix->mutable_col_blocks()) = col_blocks; - (*matrix->mutable_crsb_rows()) = crsb_rows; - (*matrix->mutable_crsb_cols()) = crsb_cols; matrix->set_storage_type(CompressedRowSparseMatrix::UNSYMMETRIC); return matrix; } diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h index 1e26f7c28..67b043e71 100644 --- a/internal/ceres/compressed_row_sparse_matrix.h +++ b/internal/ceres/compressed_row_sparse_matrix.h @@ -154,15 +154,6 @@ class CompressedRowSparseMatrix : public SparseMatrix { const std::vector& col_blocks() const { return col_blocks_; } std::vector* mutable_col_blocks() { return &col_blocks_; } - const std::vector& block_offsets() const { return block_offsets_; } - std::vector* mutable_block_offsets() { return &block_offsets_; } - - const std::vector& crsb_rows() const { return crsb_rows_; } - std::vector* mutable_crsb_rows() { return &crsb_rows_; } - - const std::vector& crsb_cols() const { return crsb_cols_; } - std::vector* mutable_crsb_cols() { return &crsb_cols_; } - // 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. @@ -217,38 +208,7 @@ class CompressedRowSparseMatrix : public SparseMatrix { static CompressedRowSparseMatrix* CreateRandomMatrix( const RandomMatrixOptions& options); - // Compute the sparsity structure of the product m.transpose() * m - // and create a CompressedRowSparseMatrix corresponding to it. - // - // Also compute a "program" vector, which for every term in the - // block outer product provides the information for the entry - // in the values array of the result matrix where it should be accumulated. - // - // This program is used by the ComputeOuterProduct function below to - // compute the outer product. - // - // Since the entries of the program are the same for rows with the - // same sparsity structure, the program only stores the result for - // one row per row block. The ComputeOuterProduct function reuses - // this information for each row in the row block. - // - // storage_type controls the form of the output matrix. It can be - // LOWER_TRIANGULAR or UPPER_TRIANGULAR. - static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram( - const CompressedRowSparseMatrix& m, - const StorageType storage_type, - std::vector* program); - - // Compute the values array for the expression m.transpose() * m, - // where the matrix used to store the result and a program have been - // created using the CreateOuterProductMatrixAndProgram function - // above. - static void ComputeOuterProduct(const CompressedRowSparseMatrix& m, - const std::vector& program, - CompressedRowSparseMatrix* result); - private: - static CompressedRowSparseMatrix* FromTripletSparseMatrix( const TripletSparseMatrix& input, bool transpose); @@ -266,18 +226,6 @@ class CompressedRowSparseMatrix : public SparseMatrix { // any way. std::vector row_blocks_; std::vector col_blocks_; - - // For outer product matrix (J' * J), we pre-compute its block - // offsets information here for fast outer product computation in - // block unit. Since the outer product matrix is symmetric, we do - // not need to distinguish row or col block. In another word, this - // is the prefix sum of row_blocks_/col_blocks_. - std::vector block_offsets_; - - // If the matrix has an underlying block structure, then it can also - // carry with it compressed row sparse block information. - std::vector crsb_rows_; - std::vector crsb_cols_; }; } // namespace internal diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc index bb8456e77..ab8f4adec 100644 --- a/internal/ceres/compressed_row_sparse_matrix_test.cc +++ b/internal/ceres/compressed_row_sparse_matrix_test.cc @@ -90,15 +90,6 @@ class CompressedRowSparseMatrixTest : public ::testing::Test { vector* col_blocks = crsm->mutable_col_blocks(); col_blocks->resize(num_cols); std::fill(col_blocks->begin(), col_blocks->end(), 1); - - // With all blocks of size 1, crsb_rows and crsb_cols are equivalent to - // rows and cols. - std::copy(crsm->rows(), - crsm->rows() + crsm->num_rows() + 1, - std::back_inserter(*crsm->mutable_crsb_rows())); - std::copy(crsm->cols(), - crsm->cols() + crsm->num_nonzeros(), - std::back_inserter(*crsm->mutable_crsb_cols())); } int num_rows; @@ -152,8 +143,6 @@ TEST_F(CompressedRowSparseMatrixTest, DeleteRows) { // Clear the row and column blocks as these are purely scalar tests. crsm->mutable_row_blocks()->clear(); crsm->mutable_col_blocks()->clear(); - crsm->mutable_crsb_rows()->clear(); - crsm->mutable_crsb_cols()->clear(); for (int i = 0; i < num_rows; ++i) { tsm->Resize(num_rows - i, num_cols); @@ -166,8 +155,6 @@ TEST_F(CompressedRowSparseMatrixTest, AppendRows) { // Clear the row and column blocks as these are purely scalar tests. crsm->mutable_row_blocks()->clear(); crsm->mutable_col_blocks()->clear(); - crsm->mutable_crsb_rows()->clear(); - crsm->mutable_crsb_cols()->clear(); for (int i = 0; i < num_rows; ++i) { TripletSparseMatrix tsm_appendage(*tsm); @@ -198,9 +185,6 @@ TEST_F(CompressedRowSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) { const vector pre_row_blocks = crsm->row_blocks(); const vector pre_col_blocks = crsm->col_blocks(); - const vector pre_crsb_rows = crsm->crsb_rows(); - const vector pre_crsb_cols = crsm->crsb_cols(); - scoped_ptr appendage( CompressedRowSparseMatrix::CreateBlockDiagonalMatrix( diagonal.get(), row_and_column_blocks)); @@ -221,22 +205,9 @@ TEST_F(CompressedRowSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) { EXPECT_EQ(expected_row_blocks, crsm->row_blocks()); EXPECT_EQ(expected_col_blocks, crsm->col_blocks()); - EXPECT_EQ(crsm->crsb_cols().size(), - pre_crsb_cols.size() + row_and_column_blocks.size()); - EXPECT_EQ(crsm->crsb_rows().size(), - pre_crsb_rows.size() + row_and_column_blocks.size()); - for (int i = 0; i < row_and_column_blocks.size(); ++i) { - EXPECT_EQ(crsm->crsb_rows()[i + pre_crsb_rows.size()], - pre_crsb_rows.back() + i + 1); - EXPECT_EQ(crsm->crsb_cols()[i + pre_crsb_cols.size()], i); - } - crsm->DeleteRows(num_diagonal_rows); EXPECT_EQ(crsm->row_blocks(), pre_row_blocks); EXPECT_EQ(crsm->col_blocks(), pre_col_blocks); - - EXPECT_EQ(crsm->crsb_rows(), pre_crsb_rows); - EXPECT_EQ(crsm->crsb_cols(), pre_crsb_cols); } TEST_F(CompressedRowSparseMatrixTest, ToDenseMatrix) { @@ -334,14 +305,6 @@ TEST(CompressedRowSparseMatrix, Transpose) { matrix.mutable_col_blocks()->push_back(4); matrix.mutable_col_blocks()->push_back(2); - matrix.mutable_crsb_rows()->push_back(0); - matrix.mutable_crsb_rows()->push_back(2); - matrix.mutable_crsb_rows()->push_back(4); - matrix.mutable_crsb_cols()->push_back(0); - matrix.mutable_crsb_cols()->push_back(1); - matrix.mutable_crsb_cols()->push_back(0); - matrix.mutable_crsb_cols()->push_back(1); - rows[0] = 0; cols[0] = 1; cols[1] = 3; @@ -392,93 +355,6 @@ TEST(CompressedRowSparseMatrix, Transpose) { EXPECT_NEAR((dense_matrix - dense_transpose.transpose()).norm(), 0.0, 1e-14); } -TEST(CompressedRowSparseMatrix, ComputeOuterProduct) { - // "Randomly generated seed." - SetRandomState(29823); - const int kMaxNumRowBlocks = 10; - const int kMaxNumColBlocks = 10; - const int kNumTrials = 10; - - // Create a random matrix, compute its outer product using Eigen and - // ComputeOuterProduct. Convert both matrices to dense matrices and - // compare their upper triangular parts. - for (int num_row_blocks = 1; num_row_blocks < kMaxNumRowBlocks; - ++num_row_blocks) { - for (int num_col_blocks = 1; num_col_blocks < kMaxNumColBlocks; - ++num_col_blocks) { - for (int trial = 0; trial < kNumTrials; ++trial) { - CompressedRowSparseMatrix::RandomMatrixOptions options; - options.num_row_blocks = num_row_blocks; - options.num_col_blocks = num_col_blocks; - options.min_row_block_size = 1; - options.max_row_block_size = 5; - options.min_col_block_size = 1; - options.max_col_block_size = 10; - options.block_density = std::max(0.1, RandDouble()); - - VLOG(2) << "num row blocks: " << options.num_row_blocks; - VLOG(2) << "num col blocks: " << options.num_col_blocks; - VLOG(2) << "min row block size: " << options.min_row_block_size; - VLOG(2) << "max row block size: " << options.max_row_block_size; - VLOG(2) << "min col block size: " << options.min_col_block_size; - VLOG(2) << "max col block size: " << options.max_col_block_size; - VLOG(2) << "block density: " << options.block_density; - - scoped_ptr random_matrix( - CompressedRowSparseMatrix::CreateRandomMatrix(options)); - - Eigen::MappedSparseMatrix mapped_random_matrix( - random_matrix->num_rows(), - random_matrix->num_cols(), - random_matrix->num_nonzeros(), - random_matrix->mutable_rows(), - random_matrix->mutable_cols(), - random_matrix->mutable_values()); - - Matrix expected_outer_product = - mapped_random_matrix.transpose() * mapped_random_matrix; - - // Use compressed row lower triangular matrix, which will then - // get mapped to a compressed column upper triangular matrix. - vector program; - scoped_ptr outer_product( - CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( - *random_matrix, - CompressedRowSparseMatrix::LOWER_TRIANGULAR, - &program)); - CompressedRowSparseMatrix::ComputeOuterProduct( - *random_matrix, program, outer_product.get()); - - EXPECT_EQ(outer_product->row_blocks(), random_matrix->col_blocks()); - EXPECT_EQ(outer_product->col_blocks(), random_matrix->col_blocks()); - - Matrix actual_outer_product = - Eigen::MappedSparseMatrix( - outer_product->num_rows(), - outer_product->num_rows(), - outer_product->num_nonzeros(), - outer_product->mutable_rows(), - outer_product->mutable_cols(), - outer_product->mutable_values()); - expected_outer_product.triangularView().setZero(); - actual_outer_product.triangularView().setZero(); - - EXPECT_EQ(actual_outer_product.rows(), actual_outer_product.cols()); - EXPECT_EQ(expected_outer_product.rows(), expected_outer_product.cols()); - EXPECT_EQ(actual_outer_product.rows(), expected_outer_product.rows()); - - const double diff_norm = - (actual_outer_product - expected_outer_product).norm() / - expected_outer_product.norm(); - EXPECT_NEAR(diff_norm, 0.0, std::numeric_limits::epsilon()) - << "expected: \n" - << expected_outer_product << "\nactual: \n" - << actual_outer_product; - } - } - } -} - TEST(CompressedRowSparseMatrix, FromTripletSparseMatrix) { TripletSparseMatrix::RandomMatrixOptions options; options.num_rows = 5;