123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181 |
- // Ceres Solver - A fast non-linear least squares minimizer
- // Copyright 2023 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/block_random_access_sparse_matrix.h"
- #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::internal {
- TEST(BlockRandomAccessSparseMatrix, GetCell) {
- ContextImpl context;
- constexpr int num_threads = 1;
- std::vector<Block> blocks;
- blocks.emplace_back(3, 0);
- blocks.emplace_back(4, 3);
- blocks.emplace_back(5, 7);
- constexpr int num_rows = 3 + 4 + 5;
- std::set<std::pair<int, int>> block_pairs;
- int num_nonzeros = 0;
- block_pairs.emplace(0, 0);
- num_nonzeros += blocks[0].size * blocks[0].size;
- block_pairs.emplace(1, 1);
- num_nonzeros += blocks[1].size * blocks[1].size;
- block_pairs.emplace(1, 2);
- num_nonzeros += blocks[1].size * blocks[2].size;
- block_pairs.emplace(0, 2);
- num_nonzeros += blocks[2].size * blocks[0].size;
- BlockRandomAccessSparseMatrix m(blocks, block_pairs, &context, num_threads);
- EXPECT_EQ(m.num_rows(), num_rows);
- EXPECT_EQ(m.num_cols(), num_rows);
- for (const auto& block_pair : block_pairs) {
- const int row_block_id = block_pair.first;
- const int col_block_id = block_pair.second;
- int row;
- int col;
- int row_stride;
- int col_stride;
- CellInfo* cell = m.GetCell(
- row_block_id, col_block_id, &row, &col, &row_stride, &col_stride);
- EXPECT_TRUE(cell != nullptr);
- EXPECT_EQ(row, 0);
- EXPECT_EQ(col, 0);
- 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].size, blocks[col_block_id].size) =
- (row_block_id + 1) * (col_block_id + 1) *
- Matrix::Ones(blocks[row_block_id].size, blocks[col_block_id].size);
- }
- const BlockSparseMatrix* bsm = m.matrix();
- EXPECT_EQ(bsm->num_nonzeros(), num_nonzeros);
- Matrix dense;
- bsm->ToDenseMatrix(&dense);
- double kTolerance = 1e-14;
- // (0, 0)
- EXPECT_NEAR(
- (dense.block(0, 0, 3, 3) - Matrix::Ones(3, 3)).norm(), 0.0, kTolerance);
- // (1, 1)
- EXPECT_NEAR((dense.block(3, 3, 4, 4) - 2 * 2 * Matrix::Ones(4, 4)).norm(),
- 0.0,
- kTolerance);
- // (1, 2)
- EXPECT_NEAR((dense.block(3, 3 + 4, 4, 5) - 2 * 3 * Matrix::Ones(4, 5)).norm(),
- 0.0,
- kTolerance);
- // (0, 2)
- EXPECT_NEAR((dense.block(0, 3 + 4, 3, 5) - 3 * 1 * Matrix::Ones(3, 5)).norm(),
- 0.0,
- kTolerance);
- // There is nothing else in the matrix besides these four blocks.
- EXPECT_NEAR(
- dense.norm(), sqrt(9. + 16. * 16. + 36. * 20. + 9. * 15.), kTolerance);
- Vector x = Vector::Ones(dense.rows());
- Vector actual_y = Vector::Zero(dense.rows());
- Vector expected_y = Vector::Zero(dense.rows());
- expected_y += dense.selfadjointView<Eigen::Upper>() * x;
- m.SymmetricRightMultiplyAndAccumulate(x.data(), actual_y.data());
- EXPECT_NEAR((expected_y - actual_y).norm(), 0.0, kTolerance)
- << "actual: " << actual_y.transpose() << "\n"
- << "expected: " << expected_y.transpose() << "matrix: \n " << dense;
- }
- // IntPairToInt64 is private, thus this fixture is needed to access and
- // test it.
- class BlockRandomAccessSparseMatrixTest : public ::testing::Test {
- public:
- void SetUp() final {
- 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, &context_, 1);
- }
- void CheckIntPairToInt64(int a, int b) {
- int64_t value = m_->IntPairToInt64(a, b);
- EXPECT_GT(value, 0) << "Overflow a = " << a << " b = " << b;
- EXPECT_GT(value, a) << "Overflow a = " << a << " b = " << b;
- EXPECT_GT(value, b) << "Overflow a = " << a << " b = " << b;
- }
- void CheckInt64ToIntPair() {
- uint64_t max_rows = m_->kRowShift;
- for (int row = max_rows - 10; row < max_rows; ++row) {
- for (int col = 0; col < 10; ++col) {
- int row_computed;
- int col_computed;
- m_->Int64ToIntPair(
- m_->IntPairToInt64(row, col), &row_computed, &col_computed);
- EXPECT_EQ(row, row_computed);
- EXPECT_EQ(col, col_computed);
- }
- }
- }
- private:
- ContextImpl context_;
- std::unique_ptr<BlockRandomAccessSparseMatrix> m_;
- };
- TEST_F(BlockRandomAccessSparseMatrixTest, IntPairToInt64Overflow) {
- CheckIntPairToInt64(std::numeric_limits<int32_t>::max(),
- std::numeric_limits<int32_t>::max());
- }
- TEST_F(BlockRandomAccessSparseMatrixTest, Int64ToIntPair) {
- CheckInt64ToIntPair();
- }
- } // namespace ceres::internal
|