diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 8a56ca01d..052e7a0a5 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -70,6 +70,7 @@ set(CERES_INTERNAL_SRC gradient_problem.cc gradient_problem_solver.cc implicit_schur_complement.cc + inner_product_computer.cc is_close.cc iterative_schur_complement_solver.cc levenberg_marquardt_strategy.cc @@ -316,6 +317,7 @@ if (BUILD_TESTING AND GFLAGS) ceres_test(graph_algorithms) ceres_test(householder_vector) ceres_test(implicit_schur_complement) + ceres_test(inner_product_computer) ceres_test(invert_psd_matrix) ceres_test(is_close) ceres_test(iterative_schur_complement_solver) diff --git a/internal/ceres/inner_product_computer.cc b/internal/ceres/inner_product_computer.cc new file mode 100644 index 000000000..052bd5acd --- /dev/null +++ b/internal/ceres/inner_product_computer.cc @@ -0,0 +1,373 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2017 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/inner_product_computer.h" + +#include + +namespace ceres { +namespace internal { +namespace { + +// Compute the product (in MATLAB notation) +// +// c(0:a_cols, 0:b_cols) = a' * b +// +// Where: +// +// a is ab_rows x a_cols +// b is ab_rows x b_cols +// c is a_cos x c_col_stride +// +// a, b and c are row-major matrices. +// +// Performance note: +// ---------------- +// +// Technically this function is a repeat of a similarly named function +// in small_blas.h but its performance is considerably better than +// that of the version there due to the way it accesses memory. +// +// TODO(sameeragarwal): Measure and tune the performance of +// small_blas.h based on the insights gained here. +EIGEN_STRONG_INLINE void MatrixTransposeMatrixMultiply(const int ab_rows, + const double* a, + const int a_cols, + const double* b, + const int b_cols, + double* c, + int c_cols) { + // Compute c as the sum of ab_rows, rank 1 outer products of the + // corresponding rows of a and b. + for (int r = 0; r < ab_rows; ++r) { + double* c_r = c; + for (int i1 = 0; i1 < a_cols; ++i1) { + const double a_v = a[i1]; + for (int i2 = 0; i2 < b_cols; ++i2) { + c_r[i2] += a_v * b[i2]; + } + c_r += c_cols; + } + a += a_cols; + b += b_cols; + } +} + +} // namespace + +// Create the CompressedRowSparseMatrix matrix that will contain the +// inner product. +// +// storage_type controls whether the result matrix contains the upper +// or the lower triangular part of the product. +// +// num_nonzeros is the number of non-zeros in the result matrix. +CompressedRowSparseMatrix* InnerProductComputer::CreateResultMatrix( + const CompressedRowSparseMatrix::StorageType storage_type, + const int num_nonzeros) { + CompressedRowSparseMatrix* matrix = + new CompressedRowSparseMatrix(m_.num_cols(), m_.num_cols(), num_nonzeros); + matrix->set_storage_type(storage_type); + + const CompressedRowBlockStructure* bs = m_.block_structure(); + const std::vector& 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; + } + + return matrix; +} + +// Given the set of product terms in the inner product, return the +// total number of non-zeros in the result and for each row block of +// the result matrix, compute the number of non-zeros in any one row +// of the row block. +int InnerProductComputer::ComputeNonzeros( + const std::vector& product_terms, + std::vector* row_nnz) { + const CompressedRowBlockStructure* bs = m_.block_structure(); + const std::vector& blocks = bs->cols; + + row_nnz->resize(blocks.size()); + std::fill(row_nnz->begin(), row_nnz->end(), 0); + + // First product term. + (*row_nnz)[product_terms[0].row] = blocks[product_terms[0].col].size; + int num_nonzeros = + blocks[product_terms[0].row].size * blocks[product_terms[0].col].size; + + // Remaining product terms. + for (int i = 1; i < product_terms.size(); ++i) { + const ProductTerm& previous = product_terms[i - 1]; + const ProductTerm& current = product_terms[i]; + + // Each (row, col) block counts only once. + // This check depends on product sorted on (row, col). + if (current.row != previous.row || current.col != previous.col) { + (*row_nnz)[current.row] += blocks[current.col].size; + num_nonzeros += blocks[current.row].size * blocks[current.col].size; + } + } + + return num_nonzeros; +} + +InnerProductComputer::InnerProductComputer(const BlockSparseMatrix& m, + const int start_row_block, + const int end_row_block) + : m_(m), start_row_block_(start_row_block), end_row_block_(end_row_block) {} + +// Compute the sparsity structure of the product m.transpose() * m +// and create a CompressedRowSparseMatrix corresponding to it. +// +// Also compute the "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. +// +// 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 Compute function reuses this information for +// each row in the row block. +// +// product_storage_type controls the form of the output matrix. It +// can be LOWER_TRIANGULAR or UPPER_TRIANGULAR. +InnerProductComputer* InnerProductComputer::Create( + const BlockSparseMatrix& m, + CompressedRowSparseMatrix::StorageType product_storage_type) { + return InnerProductComputer::Create( + m, 0, m.block_structure()->rows.size(), product_storage_type); +} + +InnerProductComputer* InnerProductComputer::Create( + const BlockSparseMatrix& m, + const int start_row_block, + const int end_row_block, + CompressedRowSparseMatrix::StorageType product_storage_type) { + CHECK(product_storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR || + product_storage_type == CompressedRowSparseMatrix::UPPER_TRIANGULAR); + CHECK_GT(m.num_nonzeros(), 0) + << "Congratulations, you found a bug in Ceres. Please report it."; + InnerProductComputer* inner_product_computer = + new InnerProductComputer(m, start_row_block, end_row_block); + inner_product_computer->Init(product_storage_type); + return inner_product_computer; +} + +void InnerProductComputer::Init( + const CompressedRowSparseMatrix::StorageType product_storage_type) { + std::vector product_terms; + const CompressedRowBlockStructure* bs = m_.block_structure(); + + // Give input matrix m in Block Sparse format + // (row_block, col_block) + // represent each block multiplication + // (row_block, col_block1)' X (row_block, col_block2) + // by its product term: + // (col_block1, col_block2, index) + for (int row_block = start_row_block_; row_block < end_row_block_; + ++row_block) { + const CompressedRow& row = bs->rows[row_block]; + for (int c1 = 0; c1 < row.cells.size(); ++c1) { + const Cell& cell1 = row.cells[c1]; + int c2_begin, c2_end; + if (product_storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR) { + c2_begin = 0; + c2_end = c1 + 1; + } else { + c2_begin = c1; + c2_end = row.cells.size(); + } + + for (int c2 = c2_begin; c2 < c2_end; ++c2) { + const Cell& cell2 = row.cells[c2]; + product_terms.push_back(InnerProductComputer::ProductTerm( + cell1.block_id, cell2.block_id, product_terms.size())); + } + } + } + + std::sort(product_terms.begin(), product_terms.end()); + ComputeOffsetsAndCreateResultMatrix(product_storage_type, product_terms); +} + +void InnerProductComputer::ComputeOffsetsAndCreateResultMatrix( + const CompressedRowSparseMatrix::StorageType product_storage_type, + const std::vector& product_terms) { + const int num_cols = m_.num_cols(); + const std::vector& col_blocks = m_.block_structure()->cols; + + std::vector row_block_nnz; + const int num_nonzeros = ComputeNonzeros(product_terms, &row_block_nnz); + + result_.reset(CreateResultMatrix(product_storage_type, num_nonzeros)); + + // Populate the row non-zero counts in the result matrix. + int* crsm_rows = result_->mutable_rows(); + crsm_rows[0] = 0; + for (int i = 0; i < col_blocks.size(); ++i) { + for (int j = 0; j < col_blocks[i].size; ++j, ++crsm_rows) { + *(crsm_rows + 1) = *crsm_rows + row_block_nnz[i]; + } + } + + // The following macro FILL_CRSM_COL_BLOCK is key to understanding + // how this class works. + // + // It does two things. + // + // Sets the value for the current term in the result_offsets_ array + // and populates the cols array of the result matrix. + // + // row_block and col_block as the names imply, refer to the row and + // column blocks of the current term. + // + // row_nnz is the number of nonzeros in the result_matrix at the + // beginning of the first row of row_block. + // + // col_nnz is the number of nonzeros in the first row of the row + // block that occur before the current column block, i.e. this is + // sum of the sizes of all the column blocks in this row block that + // came before this column block. + // + // Given these two numbers and the total number of nonzeros in this + // row (nnz_in_row), we can now populate the cols array as follows: + // + // nnz + j * nnz_in_row is the beginning of the j^th row. + // + // nnz + j * nnz_in_row + col_nnz is the beginning of the column + // block in the j^th row. + // + // nnz + j * nnz_in_row + col_nnz + k is then the j^th row and the + // k^th column of the product block, whose value is + // + // col_blocks[col_block].position + k, which is the column number of + // the k^th column of the current column block. +#define FILL_CRSM_COL_BLOCK \ + const int row_block = current->row; \ + const int col_block = current->col; \ + const int nnz_in_row = row_block_nnz[row_block]; \ + int* crsm_cols = result_->mutable_cols(); \ + result_offsets_[current->index] = nnz + col_nnz; \ + for (int j = 0; j < col_blocks[row_block].size; ++j) { \ + for (int k = 0; k < col_blocks[col_block].size; ++k) { \ + crsm_cols[nnz + j * nnz_in_row + col_nnz + k] = \ + col_blocks[col_block].position + k; \ + } \ + } + + result_offsets_.resize(product_terms.size()); + int col_nnz = 0; + int nnz = 0; + + // Process the first term. + const InnerProductComputer::ProductTerm* current = &product_terms[0]; + FILL_CRSM_COL_BLOCK; + + // Process the rest of the terms. + for (int i = 1; i < product_terms.size(); ++i) { + current = &product_terms[i]; + const InnerProductComputer::ProductTerm* previous = &product_terms[i - 1]; + + // If the current term is the same as the previous term, then it + // stores its product at the same location as the previous term. + if (previous->row == current->row && previous->col == current->col) { + result_offsets_[current->index] = result_offsets_[previous->index]; + continue; + } + + if (previous->row == current->row) { + // if the current and previous terms are in the same row block, + // then they differ in the column block, in which case advance + // col_nnz by the column size of the prevous term. + col_nnz += col_blocks[previous->col].size; + } else { + // If we have moved to a new row-block , then col_nnz is zero, + // and nnz is set to the beginning of the row block. + col_nnz = 0; + nnz += row_block_nnz[previous->row] * col_blocks[previous->row].size; + } + + FILL_CRSM_COL_BLOCK; + } +} + +// Use the results_offsets_ array to numerically compute the product +// m' * m and store it in result_. +// +// TODO(sameeragarwal): Multithreading support. +void InnerProductComputer::Compute() { + const double* m_values = m_.values(); + const CompressedRowBlockStructure* bs = m_.block_structure(); + + const CompressedRowSparseMatrix::StorageType storage_type = + result_->storage_type(); + result_->SetZero(); + double* values = result_->mutable_values(); + const int* rows = result_->rows(); + int cursor = 0; + + // Iterate row blocks. + for (int r = start_row_block_; r < end_row_block_; ++r) { + const CompressedRow& m_row = bs->rows[r]; + for (int c1 = 0; c1 < m_row.cells.size(); ++c1) { + const Cell& cell1 = m_row.cells[c1]; + const int c1_size = bs->cols[cell1.block_id].size; + const int row_nnz = rows[bs->cols[cell1.block_id].position + 1] - + rows[bs->cols[cell1.block_id].position]; + + int c2_begin, c2_end; + if (storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR) { + c2_begin = 0; + c2_end = c1 + 1; + } else { + c2_begin = c1; + c2_end = m_row.cells.size(); + } + + for (int c2 = c2_begin; c2 < c2_end; ++c2, ++cursor) { + const Cell& cell2 = m_row.cells[c2]; + const int c2_size = bs->cols[cell2.block_id].size; + MatrixTransposeMatrixMultiply(m_row.block.size, + m_values + cell1.position, c1_size, + m_values + cell2.position, c2_size, + values + result_offsets_[cursor], + row_nnz); + } + } + } + + CHECK_EQ(cursor, result_offsets_.size()); +} + +} // namespace internal +} // namespace ceres diff --git a/internal/ceres/inner_product_computer.h b/internal/ceres/inner_product_computer.h new file mode 100644 index 000000000..d206707d3 --- /dev/null +++ b/internal/ceres/inner_product_computer.h @@ -0,0 +1,157 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2017 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#ifndef CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_ +#define CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_ + +#include + +#include "ceres/block_sparse_matrix.h" +#include "ceres/compressed_row_sparse_matrix.h" +#include "ceres/internal/scoped_ptr.h" + +namespace ceres { +namespace internal { + +// This class is used to repeatedly compute the inner product +// +// result = m' * m +// +// where the sparsity structure of m remains constant across calls. +// +// Upon creation, the class computes and caches information needed to +// compute v, and then uses it to efficiently compute the product +// every time InnerProductComputer::Compute is called. +// +// See sparse_normal_cholesky_solver.cc for example usage. +// +// Note that the result matrix is a block upper or lower-triangular +// matrix, i.e., it will contain entries in the upper or lower +// triangular part of the matrix corresponding to the block that occur +// along its diagonal. +// +// This is not a problem as sparse linear algebra libraries can ignore +// these entries with ease and the space used is minimal/linear in the +// size of the matrices. +class InnerProductComputer { + public: + // Factory + // + // m is the input matrix + // + // Since m' * m is a symmetric matrix, we only compute half of the + // matrix and the value of storage_type which must be + // UPPER_TRIANGULAR or LOWER_TRIANGULAR determines which half is + // computed. + // + // The user must ensure that the matrix m is valid for the life time + // of this object. + static InnerProductComputer* Create( + const BlockSparseMatrix& m, + CompressedRowSparseMatrix::StorageType storage_type); + + // This factory method allows the user control over range of row + // blocks of m that should be used to compute the inner product. + // + // a = m(start_row_block : end_row_block, :); + // result = a' * a; + static InnerProductComputer* Create( + const BlockSparseMatrix& m, + int start_row_block, + int end_row_block, + CompressedRowSparseMatrix::StorageType storage_type); + + // Update result_ to be numerically equal to m' * m. + void Compute(); + + // Accessors for the result containing the inner product. + // + // Compute must be called before accessing this result for + // the first time. + const CompressedRowSparseMatrix& result() const { return *result_; } + CompressedRowSparseMatrix* mutable_result() const { return result_.get(); } + + private: + // A ProductTerm is a term in the block inner 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; + }; + + InnerProductComputer(const BlockSparseMatrix& m, + int start_row_block, + int end_row_block); + + void Init(CompressedRowSparseMatrix::StorageType storage_type); + + CompressedRowSparseMatrix* CreateResultMatrix( + const CompressedRowSparseMatrix::StorageType storage_type, + int num_nonzeros); + + int ComputeNonzeros(const std::vector& product_terms, + std::vector* row_block_nnz); + + void ComputeOffsetsAndCreateResultMatrix( + const CompressedRowSparseMatrix::StorageType storage_type, + const std::vector& product_terms); + + const BlockSparseMatrix& m_; + const int start_row_block_; + const int end_row_block_; + scoped_ptr result_; + + // For each term in the inner product, result_offsets_ contains the + // location in the values array of the result_ matrix where it + // should be stored. + // + // This is the principal look up table that allows this class to + // compute the inner product fast. + std::vector result_offsets_; +}; + +} // namespace internal +} // namespace ceres + +#endif // CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_ diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc new file mode 100644 index 000000000..a4c81bcdc --- /dev/null +++ b/internal/ceres/inner_product_computer_test.cc @@ -0,0 +1,236 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2017 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/inner_product_computer.h" + +#include +#include "ceres/block_sparse_matrix.h" +#include "ceres/internal/eigen.h" +#include "ceres/internal/scoped_ptr.h" +#include "ceres/random.h" +#include "ceres/triplet_sparse_matrix.h" +#include "glog/logging.h" +#include "gtest/gtest.h" + +#include "Eigen/SparseCore" + +namespace ceres { +namespace internal { + +template +void CompareTriangularPartOfMatrices(const Matrix& expected, + const Matrix& actual) { + EXPECT_EQ(actual.rows(), actual.cols()); + EXPECT_EQ(expected.rows(), expected.cols()); + EXPECT_EQ(actual.rows(), expected.rows()); + + const Matrix expected_t = expected.triangularView(); + const Matrix actual_t = actual.triangularView(); + + // TODO(sameeragarwal): Foo + CHECK_LE((expected_t - actual_t).norm() / actual_t.norm(), + 100 * std::numeric_limits::epsilon()) + << "expected: \n" + << expected_t << "\nactual: \n" + << actual_t; +} + +#define COMPUTE_AND_COMPARE \ + { \ + inner_product_computer->Compute(); \ + CompressedRowSparseMatrix* actual_product_crsm = \ + inner_product_computer->mutable_result(); \ + Matrix actual_inner_product = \ + Eigen::MappedSparseMatrix( \ + actual_product_crsm->num_rows(), \ + actual_product_crsm->num_rows(), \ + actual_product_crsm->num_nonzeros(), \ + actual_product_crsm->mutable_rows(), \ + actual_product_crsm->mutable_cols(), \ + actual_product_crsm->mutable_values()); \ + (actual_product_crsm->storage_type() == \ + CompressedRowSparseMatrix::LOWER_TRIANGULAR) \ + ? CompareTriangularPartOfMatrices( \ + expected_inner_product, actual_inner_product) \ + : CompareTriangularPartOfMatrices( \ + expected_inner_product, actual_inner_product); \ + } + + +TEST(InnerProductComputer, NormalOperation) { + // "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) { + BlockSparseMatrix::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( + BlockSparseMatrix::CreateRandomMatrix(options)); + + TripletSparseMatrix tsm(random_matrix->num_rows(), + random_matrix->num_cols(), + random_matrix->num_nonzeros()); + random_matrix->ToTripletSparseMatrix(&tsm); + std::vector > triplets; + for (int i = 0; i < tsm.num_nonzeros(); ++i) { + triplets.push_back(Eigen::Triplet( + tsm.rows()[i], tsm.cols()[i], tsm.values()[i])); + } + Eigen::SparseMatrix eigen_random_matrix( + random_matrix->num_rows(), random_matrix->num_cols()); + eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end()); + Matrix expected_inner_product = + eigen_random_matrix.transpose() * eigen_random_matrix; + + scoped_ptr inner_product_computer; + + inner_product_computer.reset(InnerProductComputer::Create( + *random_matrix, CompressedRowSparseMatrix::LOWER_TRIANGULAR)); + COMPUTE_AND_COMPARE; + inner_product_computer.reset(InnerProductComputer::Create( + *random_matrix, CompressedRowSparseMatrix::UPPER_TRIANGULAR)); + COMPUTE_AND_COMPARE; + + } + } + } +} + + +TEST(InnerProductComputer, SubMatrix) { + // "Randomly generated seed." + SetRandomState(29823); + const int kNumRowBlocks = 10; + const int kNumColBlocks = 20; + const int kNumTrials = 5; + + // Create a random matrix, compute its outer product using Eigen and + // ComputeInnerProductComputer. Convert both matrices to dense matrices and + // compare their upper triangular parts. + for (int trial = 0; trial < kNumTrials; ++trial) { + BlockSparseMatrix::RandomMatrixOptions options; + options.num_row_blocks = kNumRowBlocks; + options.num_col_blocks = kNumColBlocks; + 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( + BlockSparseMatrix::CreateRandomMatrix(options)); + + const std::vector& row_blocks = + random_matrix->block_structure()->rows; + const int num_row_blocks = row_blocks.size(); + + for (int start_row_block = 0; start_row_block < num_row_blocks - 1; + ++start_row_block) { + for (int end_row_block = start_row_block + 1; + end_row_block < num_row_blocks; + ++end_row_block) { + const int start_row = row_blocks[start_row_block].block.position; + const int end_row = row_blocks[end_row_block].block.position; + + TripletSparseMatrix tsm(random_matrix->num_rows(), + random_matrix->num_cols(), + random_matrix->num_nonzeros()); + random_matrix->ToTripletSparseMatrix(&tsm); + std::vector > triplets; + for (int i = 0; i < tsm.num_nonzeros(); ++i) { + if (tsm.rows()[i] >= start_row && tsm.rows()[i] < end_row) { + triplets.push_back(Eigen::Triplet( + tsm.rows()[i], tsm.cols()[i], tsm.values()[i])); + } + } + + Eigen::SparseMatrix eigen_random_matrix( + random_matrix->num_rows(), random_matrix->num_cols()); + eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end()); + + Matrix expected_inner_product = + eigen_random_matrix.transpose() * eigen_random_matrix; + + scoped_ptr inner_product_computer; + inner_product_computer.reset(InnerProductComputer::Create( + *random_matrix, + start_row_block, + end_row_block, + CompressedRowSparseMatrix::LOWER_TRIANGULAR)); + COMPUTE_AND_COMPARE; + inner_product_computer.reset(InnerProductComputer::Create( + *random_matrix, + start_row_block, + end_row_block, + CompressedRowSparseMatrix::UPPER_TRIANGULAR)); + COMPUTE_AND_COMPARE; + + } + } + } +} + +#undef COMPUTE_AND_COMPARE + +} // namespace internal +} // namespace ceres diff --git a/jni/Android.mk b/jni/Android.mk index 20b484d46..551cd45b0 100644 --- a/jni/Android.mk +++ b/jni/Android.mk @@ -154,6 +154,7 @@ LOCAL_SRC_FILES := $(CERES_SRC_PATH)/array_utils.cc \ $(CERES_SRC_PATH)/gradient_problem_solver.cc \ $(CERES_SRC_PATH)/is_close.cc \ $(CERES_SRC_PATH)/implicit_schur_complement.cc \ + $(CERES_SRC_PATH)/inner_product_computer.cc \ $(CERES_SRC_PATH)/iterative_schur_complement_solver.cc \ $(CERES_SRC_PATH)/lapack.cc \ $(CERES_SRC_PATH)/levenberg_marquardt_strategy.cc \