From 04899645cc75160ed74fe88edbb7c9f516d21360 Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Wed, 10 Aug 2022 09:55:43 -0700 Subject: [PATCH] LinearOperator::FooMultiply -> LinearOperator::FooMultiplyAndAccumulate These methods were historically poorly named and every time I read code I get confused whether they are just multiplying or multiplying and adding. Clarifying them also gives us the changce to introduce RightMultiply and LeftMultiply methods in the base class which will simplify a number call sites in a subsequent CL. Fixes https://github.com/ceres-solver/ceres-solver/issues/855 Change-Id: Ice4fb483f1acd02527a6dd753ef0c5a66037f4b0 --- internal/ceres/block_jacobi_preconditioner.cc | 6 ++-- internal/ceres/block_jacobi_preconditioner.h | 2 +- .../block_random_access_diagonal_matrix.cc | 4 +-- .../block_random_access_diagonal_matrix.h | 2 +- ...lock_random_access_diagonal_matrix_test.cc | 4 +-- .../block_random_access_sparse_matrix.cc | 4 +-- .../ceres/block_random_access_sparse_matrix.h | 2 +- .../block_random_access_sparse_matrix_test.cc | 2 +- internal/ceres/block_sparse_matrix.cc | 6 ++-- internal/ceres/block_sparse_matrix.h | 4 +-- internal/ceres/block_sparse_matrix_test.cc | 26 +++++++------- internal/ceres/cgnr_solver.cc | 8 ++--- .../ceres/compressed_row_sparse_matrix.cc | 16 +++++---- internal/ceres/compressed_row_sparse_matrix.h | 4 +-- .../compressed_row_sparse_matrix_test.cc | 25 ++++++------- internal/ceres/conjugate_gradients_solver.h | 15 ++++---- internal/ceres/context_impl.h | 2 +- internal/ceres/cuda_buffer.h | 5 +-- internal/ceres/cuda_kernels_test.cc | 28 +++++---------- internal/ceres/dense_sparse_matrix.cc | 6 ++-- internal/ceres/dense_sparse_matrix.h | 4 +-- internal/ceres/dense_sparse_matrix_test.cc | 20 +++++------ internal/ceres/dogleg_strategy.cc | 8 ++--- .../dynamic_sparse_normal_cholesky_solver.cc | 2 +- ...amic_sparse_normal_cholesky_solver_test.cc | 2 +- internal/ceres/implicit_schur_complement.cc | 35 ++++++++++--------- internal/ceres/implicit_schur_complement.h | 14 ++++---- .../ceres/implicit_schur_complement_test.cc | 2 +- internal/ceres/iterative_refiner.cc | 2 +- internal/ceres/iterative_refiner_test.cc | 6 ++-- internal/ceres/line_search_direction.cc | 4 +-- internal/ceres/linear_operator.h | 12 +++---- internal/ceres/low_rank_inverse_hessian.cc | 4 +-- internal/ceres/low_rank_inverse_hessian.h | 8 ++--- internal/ceres/partitioned_matrix_view.h | 18 +++++----- internal/ceres/partitioned_matrix_view_impl.h | 16 ++++----- .../ceres/partitioned_matrix_view_test.cc | 20 +++++------ internal/ceres/preconditioner.cc | 6 ++-- internal/ceres/preconditioner.h | 17 ++++----- internal/ceres/schur_complement_solver.cc | 8 ++--- internal/ceres/schur_jacobi_preconditioner.cc | 6 ++-- internal/ceres/schur_jacobi_preconditioner.h | 4 +-- internal/ceres/sparse_matrix.h | 5 +-- .../ceres/sparse_normal_cholesky_solver.cc | 2 +- .../sparse_normal_cholesky_solver_test.cc | 2 +- internal/ceres/subset_preconditioner.cc | 3 +- internal/ceres/subset_preconditioner.h | 2 +- internal/ceres/subset_preconditioner_test.cc | 2 +- internal/ceres/triplet_sparse_matrix.cc | 6 ++-- internal/ceres/triplet_sparse_matrix.h | 5 +-- internal/ceres/triplet_sparse_matrix_test.cc | 2 +- internal/ceres/trust_region_minimizer.cc | 3 +- .../ceres/visibility_based_preconditioner.cc | 4 +-- .../ceres/visibility_based_preconditioner.h | 4 +-- .../visibility_based_preconditioner_test.cc | 2 +- 55 files changed, 221 insertions(+), 210 deletions(-) diff --git a/internal/ceres/block_jacobi_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc index 07d6cd38a..fdba2b857 100644 --- a/internal/ceres/block_jacobi_preconditioner.cc +++ b/internal/ceres/block_jacobi_preconditioner.cc @@ -90,9 +90,9 @@ bool BlockJacobiPreconditioner::UpdateImpl(const BlockSparseMatrix& A, return true; } -void BlockJacobiPreconditioner::RightMultiply(const double* x, - double* y) const { - m_->RightMultiply(x, y); +void BlockJacobiPreconditioner::RightMultiplyAndAccumulate(const double* x, + double* y) const { + m_->RightMultiplyAndAccumulate(x, y); } } // namespace ceres::internal diff --git a/internal/ceres/block_jacobi_preconditioner.h b/internal/ceres/block_jacobi_preconditioner.h index 0928919f1..7728eb931 100644 --- a/internal/ceres/block_jacobi_preconditioner.h +++ b/internal/ceres/block_jacobi_preconditioner.h @@ -64,7 +64,7 @@ class CERES_NO_EXPORT BlockJacobiPreconditioner ~BlockJacobiPreconditioner() override; // Preconditioner interface - void RightMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; int num_rows() const final { return m_->num_rows(); } int num_cols() const final { return m_->num_rows(); } const BlockRandomAccessDiagonalMatrix& matrix() const { return *m_; } diff --git a/internal/ceres/block_random_access_diagonal_matrix.cc b/internal/ceres/block_random_access_diagonal_matrix.cc index 5a9f772ae..006713fb5 100644 --- a/internal/ceres/block_random_access_diagonal_matrix.cc +++ b/internal/ceres/block_random_access_diagonal_matrix.cc @@ -130,8 +130,8 @@ void BlockRandomAccessDiagonalMatrix::Invert() { } } -void BlockRandomAccessDiagonalMatrix::RightMultiply(const double* x, - double* y) const { +void BlockRandomAccessDiagonalMatrix::RightMultiplyAndAccumulate( + const double* x, double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); const double* values = tsm_->values(); diff --git a/internal/ceres/block_random_access_diagonal_matrix.h b/internal/ceres/block_random_access_diagonal_matrix.h index 45e3d02d7..2a726a0ed 100644 --- a/internal/ceres/block_random_access_diagonal_matrix.h +++ b/internal/ceres/block_random_access_diagonal_matrix.h @@ -75,7 +75,7 @@ class CERES_NO_EXPORT BlockRandomAccessDiagonalMatrix void Invert(); // y += S * x - void RightMultiply(const double* x, double* y) const; + void RightMultiplyAndAccumulate(const double* x, double* y) const; // Since the matrix is square, num_rows() == num_cols(). int num_rows() const final { return tsm_->num_rows(); } diff --git a/internal/ceres/block_random_access_diagonal_matrix_test.cc b/internal/ceres/block_random_access_diagonal_matrix_test.cc index 42a309f46..37e1f8833 100644 --- a/internal/ceres/block_random_access_diagonal_matrix_test.cc +++ b/internal/ceres/block_random_access_diagonal_matrix_test.cc @@ -134,7 +134,7 @@ TEST_F(BlockRandomAccessDiagonalMatrixTest, MatrixContents) { kTolerance); } -TEST_F(BlockRandomAccessDiagonalMatrixTest, RightMultiply) { +TEST_F(BlockRandomAccessDiagonalMatrixTest, RightMultiplyAndAccumulate) { double kTolerance = 1e-14; const TripletSparseMatrix* tsm = m_->matrix(); Matrix dense; @@ -142,7 +142,7 @@ TEST_F(BlockRandomAccessDiagonalMatrixTest, RightMultiply) { Vector x = Vector::Random(dense.rows()); Vector expected_y = dense * x; Vector actual_y = Vector::Zero(dense.rows()); - m_->RightMultiply(x.data(), actual_y.data()); + m_->RightMultiplyAndAccumulate(x.data(), actual_y.data()); EXPECT_NEAR((expected_y - actual_y).norm(), 0, kTolerance); } diff --git a/internal/ceres/block_random_access_sparse_matrix.cc b/internal/ceres/block_random_access_sparse_matrix.cc index 3ae0cbaf4..2df4c71f2 100644 --- a/internal/ceres/block_random_access_sparse_matrix.cc +++ b/internal/ceres/block_random_access_sparse_matrix.cc @@ -147,8 +147,8 @@ void BlockRandomAccessSparseMatrix::SetZero() { } } -void BlockRandomAccessSparseMatrix::SymmetricRightMultiply(const double* x, - double* y) const { +void BlockRandomAccessSparseMatrix::SymmetricRightMultiplyAndAccumulate( + const double* x, double* y) const { for (const auto& cell_position_and_data : cell_values_) { const int row = cell_position_and_data.first.first; const int row_block_size = blocks_[row]; diff --git a/internal/ceres/block_random_access_sparse_matrix.h b/internal/ceres/block_random_access_sparse_matrix.h index 882292c53..fe2b13c88 100644 --- a/internal/ceres/block_random_access_sparse_matrix.h +++ b/internal/ceres/block_random_access_sparse_matrix.h @@ -83,7 +83,7 @@ class CERES_NO_EXPORT BlockRandomAccessSparseMatrix // matrix is stored. // // y += S * x - void SymmetricRightMultiply(const double* x, double* y) const; + void SymmetricRightMultiplyAndAccumulate(const double* x, double* y) const; // Since the matrix is square, num_rows() == num_cols(). int num_rows() const final { return tsm_->num_rows(); } diff --git a/internal/ceres/block_random_access_sparse_matrix_test.cc b/internal/ceres/block_random_access_sparse_matrix_test.cc index 7224b65f9..605bfbeb6 100644 --- a/internal/ceres/block_random_access_sparse_matrix_test.cc +++ b/internal/ceres/block_random_access_sparse_matrix_test.cc @@ -127,7 +127,7 @@ TEST(BlockRandomAccessSparseMatrix, GetCell) { Vector expected_y = Vector::Zero(dense.rows()); expected_y += dense.selfadjointView() * x; - m.SymmetricRightMultiply(x.data(), actual_y.data()); + 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; diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index 1bfa343ec..ae6bd3a75 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc @@ -88,7 +88,8 @@ void BlockSparseMatrix::SetZero() { std::fill(values_.get(), values_.get() + num_nonzeros_, 0.0); } -void BlockSparseMatrix::RightMultiply(const double* x, double* y) const { +void BlockSparseMatrix::RightMultiplyAndAccumulate(const double* x, + double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); @@ -110,7 +111,8 @@ void BlockSparseMatrix::RightMultiply(const double* x, double* y) const { } } -void BlockSparseMatrix::LeftMultiply(const double* x, double* y) const { +void BlockSparseMatrix::LeftMultiplyAndAccumulate(const double* x, + double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h index da6b641f7..7cef18deb 100644 --- a/internal/ceres/block_sparse_matrix.h +++ b/internal/ceres/block_sparse_matrix.h @@ -71,8 +71,8 @@ class CERES_NO_EXPORT BlockSparseMatrix final : public SparseMatrix { // Implementation of SparseMatrix interface. void SetZero() final; - void RightMultiply(const double* x, double* y) const final; - void LeftMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; + void LeftMultiplyAndAccumulate(const double* x, double* y) const final; void SquaredColumnNorm(double* x) const final; void ScaleColumns(const double* scale) final; void ToCRSMatrix(CRSMatrix* matrix) const; diff --git a/internal/ceres/block_sparse_matrix_test.cc b/internal/ceres/block_sparse_matrix_test.cc index 7fab13ae7..4b02abf82 100644 --- a/internal/ceres/block_sparse_matrix_test.cc +++ b/internal/ceres/block_sparse_matrix_test.cc @@ -148,26 +148,26 @@ TEST_F(BlockSparseMatrixTest, SetZeroTest) { EXPECT_EQ(13, A_->num_nonzeros()); } -TEST_F(BlockSparseMatrixTest, RightMultiplyTest) { +TEST_F(BlockSparseMatrixTest, RightMultiplyAndAccumulateTest) { Vector y_a = Vector::Zero(A_->num_rows()); Vector y_b = Vector::Zero(A_->num_rows()); for (int i = 0; i < A_->num_cols(); ++i) { Vector x = Vector::Zero(A_->num_cols()); x[i] = 1.0; - A_->RightMultiply(x.data(), y_a.data()); - B_->RightMultiply(x.data(), y_b.data()); + A_->RightMultiplyAndAccumulate(x.data(), y_a.data()); + B_->RightMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_LT((y_a - y_b).norm(), 1e-12); } } -TEST_F(BlockSparseMatrixTest, LeftMultiplyTest) { +TEST_F(BlockSparseMatrixTest, LeftMultiplyAndAccumulateTest) { Vector y_a = Vector::Zero(A_->num_cols()); Vector y_b = Vector::Zero(A_->num_cols()); for (int i = 0; i < A_->num_rows(); ++i) { Vector x = Vector::Zero(A_->num_rows()); x[i] = 1.0; - A_->LeftMultiply(x.data(), y_a.data()); - B_->LeftMultiply(x.data(), y_b.data()); + A_->LeftMultiplyAndAccumulate(x.data(), y_a.data()); + B_->LeftMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_LT((y_a - y_b).norm(), 1e-12); } } @@ -210,8 +210,8 @@ TEST_F(BlockSparseMatrixTest, AppendRows) { y_a.setZero(); y_b.setZero(); - A_->RightMultiply(x.data(), y_a.data()); - B_->RightMultiply(x.data(), y_b.data()); + A_->RightMultiplyAndAccumulate(x.data(), y_a.data()); + B_->RightMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_LT((y_a - y_b).norm(), 1e-12); } } @@ -237,8 +237,8 @@ TEST_F(BlockSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) { y_a.setZero(); y_b.setZero(); - A_->RightMultiply(x.data(), y_a.data()); - B_->RightMultiply(x.data(), y_b.data()); + A_->RightMultiplyAndAccumulate(x.data(), y_a.data()); + B_->RightMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_LT((y_a.head(B_->num_rows()) - y_b.head(B_->num_rows())).norm(), 1e-12); Vector expected_tail = Vector::Zero(A_->num_cols()); @@ -258,8 +258,8 @@ TEST_F(BlockSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) { y_a.setZero(); y_b.setZero(); - A_->RightMultiply(x.data(), y_a.data()); - B_->RightMultiply(x.data(), y_b.data()); + A_->RightMultiplyAndAccumulate(x.data(), y_a.data()); + B_->RightMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_LT((y_a - y_b).norm(), 1e-12); } } @@ -287,7 +287,7 @@ TEST(BlockSparseMatrix, CreateDiagonalMatrix) { EXPECT_EQ(m->num_rows(), m->num_cols()); Vector x = Vector::Ones(num_cols); Vector y = Vector::Zero(num_cols); - m->RightMultiply(x.data(), y.data()); + m->RightMultiplyAndAccumulate(x.data(), y.data()); for (int i = 0; i < num_cols; ++i) { EXPECT_NEAR(y[i], diagonal[i], std::numeric_limits::epsilon()); } diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc index f79b897e3..99d53755b 100644 --- a/internal/ceres/cgnr_solver.cc +++ b/internal/ceres/cgnr_solver.cc @@ -85,12 +85,12 @@ class CERES_NO_EXPORT CgnrLinearOperator final CgnrLinearOperator(const LinearOperator& A, const double* D) : A_(A), D_(D), z_(Vector::Zero(A.num_rows())) {} - void RightMultiply(const Vector& x, Vector& y) final { + void RightMultiplyAndAccumulate(const Vector& x, Vector& y) final { // z = Ax // y = y + Atz z_.setZero(); - A_.RightMultiply(x, z_); - A_.LeftMultiply(z_, y); + A_.RightMultiplyAndAccumulate(x, z_); + A_.LeftMultiplyAndAccumulate(z_, y); // y = y + DtDx if (D_ != nullptr) { @@ -159,7 +159,7 @@ LinearSolver::Summary CgnrSolver::SolveImpl( // rhs = Atb. Vector rhs(A->num_cols()); rhs.setZero(); - A->LeftMultiply(b, rhs.data()); + A->LeftMultiplyAndAccumulate(b, rhs.data()); cg_solution_ = Vector::Zero(A->num_cols()); for (int i = 0; i < 4; ++i) { diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc index 9c7856599..574b0c022 100644 --- a/internal/ceres/compressed_row_sparse_matrix.cc +++ b/internal/ceres/compressed_row_sparse_matrix.cc @@ -278,10 +278,10 @@ void CompressedRowSparseMatrix::SetZero() { std::fill(values_.begin(), values_.end(), 0); } -// TODO(sameeragarwal): Make RightMultiply and LeftMultiply -// block-aware for higher performance. -void CompressedRowSparseMatrix::RightMultiply(const double* x, - double* y) const { +// TODO(sameeragarwal): Make RightMultiplyAndAccumulate and +// LeftMultiplyAndAccumulate block-aware for higher performance. +void CompressedRowSparseMatrix::RightMultiplyAndAccumulate(const double* x, + double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); @@ -342,7 +342,8 @@ void CompressedRowSparseMatrix::RightMultiply(const double* x, } } -void CompressedRowSparseMatrix::LeftMultiply(const double* x, double* y) const { +void CompressedRowSparseMatrix::LeftMultiplyAndAccumulate(const double* x, + double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); @@ -353,8 +354,9 @@ void CompressedRowSparseMatrix::LeftMultiply(const double* x, double* y) const { } } } else { - // Since the matrix is symmetric, LeftMultiply = RightMultiply. - RightMultiply(x, y); + // Since the matrix is symmetric, LeftMultiplyAndAccumulate = + // RightMultiplyAndAccumulate. + RightMultiplyAndAccumulate(x, y); } } diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h index 1d1ac956a..25800455a 100644 --- a/internal/ceres/compressed_row_sparse_matrix.h +++ b/internal/ceres/compressed_row_sparse_matrix.h @@ -100,8 +100,8 @@ class CERES_NO_EXPORT CompressedRowSparseMatrix : public SparseMatrix { // SparseMatrix interface. ~CompressedRowSparseMatrix() override; void SetZero() final; - void RightMultiply(const double* x, double* y) const final; - void LeftMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; + void LeftMultiplyAndAccumulate(const double* x, double* y) const final; void SquaredColumnNorm(double* x) const final; void ScaleColumns(const double* scale) final; void ToDenseMatrix(Matrix* dense_matrix) const final; diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc index 0f1b94890..42f5498d0 100644 --- a/internal/ceres/compressed_row_sparse_matrix_test.cc +++ b/internal/ceres/compressed_row_sparse_matrix_test.cc @@ -64,8 +64,8 @@ static void CompareMatrices(const SparseMatrix* a, const SparseMatrix* b) { Vector y_a = Vector::Zero(num_rows); Vector y_b = Vector::Zero(num_rows); - a->RightMultiply(x.data(), y_a.data()); - b->RightMultiply(x.data(), y_b.data()); + a->RightMultiplyAndAccumulate(x.data(), y_a.data()); + b->RightMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_EQ((y_a - y_b).norm(), 0); } } @@ -237,13 +237,13 @@ TEST(CompressedRowSparseMatrix, CreateBlockDiagonalMatrix) { x.setOnes(); y.setZero(); - matrix->RightMultiply(x.data(), y.data()); + matrix->RightMultiplyAndAccumulate(x.data(), y.data()); for (int i = 0; i < diagonal.size(); ++i) { EXPECT_EQ(y[i], diagonal[i]); } y.setZero(); - matrix->LeftMultiply(x.data(), y.data()); + matrix->LeftMultiplyAndAccumulate(x.data(), y.data()); for (int i = 0; i < diagonal.size(); ++i) { EXPECT_EQ(y[i], diagonal[i]); } @@ -396,9 +396,10 @@ static std::string ParamInfoToString(testing::TestParamInfo info) { return "UNSYMMETRIC"; } -class RightMultiplyTest : public ::testing::TestWithParam {}; +class RightMultiplyAndAccumulateTest : public ::testing::TestWithParam { +}; -TEST_P(RightMultiplyTest, _) { +TEST_P(RightMultiplyAndAccumulateTest, _) { const int kMinNumBlocks = 1; const int kMaxNumBlocks = 10; const int kMinBlockSize = 1; @@ -429,7 +430,7 @@ TEST_P(RightMultiplyTest, _) { Vector actual_y(num_rows); actual_y.setZero(); - matrix->RightMultiply(x.data(), actual_y.data()); + matrix->RightMultiplyAndAccumulate(x.data(), actual_y.data()); Matrix dense; matrix->ToDenseMatrix(&dense); @@ -460,15 +461,15 @@ TEST_P(RightMultiplyTest, _) { INSTANTIATE_TEST_SUITE_P( CompressedRowSparseMatrix, - RightMultiplyTest, + RightMultiplyAndAccumulateTest, ::testing::Values(CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR, CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR, CompressedRowSparseMatrix::StorageType::UNSYMMETRIC), ParamInfoToString); -class LeftMultiplyTest : public ::testing::TestWithParam {}; +class LeftMultiplyAndAccumulateTest : public ::testing::TestWithParam {}; -TEST_P(LeftMultiplyTest, _) { +TEST_P(LeftMultiplyAndAccumulateTest, _) { const int kMinNumBlocks = 1; const int kMaxNumBlocks = 10; const int kMinBlockSize = 1; @@ -499,7 +500,7 @@ TEST_P(LeftMultiplyTest, _) { Vector actual_y(num_cols); actual_y.setZero(); - matrix->LeftMultiply(x.data(), actual_y.data()); + matrix->LeftMultiplyAndAccumulate(x.data(), actual_y.data()); Matrix dense; matrix->ToDenseMatrix(&dense); @@ -530,7 +531,7 @@ TEST_P(LeftMultiplyTest, _) { INSTANTIATE_TEST_SUITE_P( CompressedRowSparseMatrix, - LeftMultiplyTest, + LeftMultiplyAndAccumulateTest, ::testing::Values(CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR, CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR, CompressedRowSparseMatrix::StorageType::UNSYMMETRIC), diff --git a/internal/ceres/conjugate_gradients_solver.h b/internal/ceres/conjugate_gradients_solver.h index 6254d2cb3..93f9e25e9 100644 --- a/internal/ceres/conjugate_gradients_solver.h +++ b/internal/ceres/conjugate_gradients_solver.h @@ -55,7 +55,8 @@ template class ConjugateGradientsLinearOperator { public: ~ConjugateGradientsLinearOperator() = default; - virtual void RightMultiply(const DenseVectorType& x, DenseVectorType& y) = 0; + virtual void RightMultiplyAndAccumulate(const DenseVectorType& x, + DenseVectorType& y) = 0; }; // Adapter class that makes LinearOperator appear like an instance of @@ -65,8 +66,8 @@ class LinearOperatorAdapter : public ConjugateGradientsLinearOperator { LinearOperatorAdapter(LinearOperator& linear_operator) : linear_operator_(linear_operator) {} - void RightMultiply(const Vector& x, Vector& y) final { - linear_operator_.RightMultiply(x, y); + void RightMultiplyAndAccumulate(const Vector& x, Vector& y) final { + linear_operator_.RightMultiplyAndAccumulate(x, y); } private: @@ -135,7 +136,7 @@ LinearSolver::Summary ConjugateGradientsSolver( const double tol_r = options.r_tolerance * norm_b; SetZero(tmp); - lhs.RightMultiply(solution, tmp); + lhs.RightMultiplyAndAccumulate(solution, tmp); // r = rhs - tmp Axpby(1.0, rhs, -1.0, tmp, r); @@ -157,7 +158,7 @@ LinearSolver::Summary ConjugateGradientsSolver( for (summary.num_iterations = 1;; ++summary.num_iterations) { SetZero(z); - preconditioner.RightMultiply(r, z); + preconditioner.RightMultiplyAndAccumulate(r, z); double last_rho = rho; // rho = r.dot(z); @@ -188,7 +189,7 @@ LinearSolver::Summary ConjugateGradientsSolver( DenseVectorType& q = z; SetZero(q); - lhs.RightMultiply(p, q); + lhs.RightMultiplyAndAccumulate(p, q); const double pq = Dot(p, q); if ((pq <= 0) || std::isinf(pq)) { summary.termination_type = LinearSolverTerminationType::NO_CONVERGENCE; @@ -224,7 +225,7 @@ LinearSolver::Summary ConjugateGradientsSolver( // double the complexity of the CG algorithm. if (summary.num_iterations % options.residual_reset_period == 0) { SetZero(tmp); - lhs.RightMultiply(solution, tmp); + lhs.RightMultiplyAndAccumulate(solution, tmp); Axpby(1.0, rhs, -1.0, tmp, r); // r = rhs - tmp; } else { diff --git a/internal/ceres/context_impl.h b/internal/ceres/context_impl.h index 3324b52c9..d4bd436ac 100644 --- a/internal/ceres/context_impl.h +++ b/internal/ceres/context_impl.h @@ -45,8 +45,8 @@ #ifndef CERES_NO_CUDA #include "cublas_v2.h" #include "cuda_runtime.h" -#include "cusparse.h" #include "cusolverDn.h" +#include "cusparse.h" #endif // CERES_NO_CUDA #ifdef CERES_USE_CXX_THREADS diff --git a/internal/ceres/cuda_buffer.h b/internal/ceres/cuda_buffer.h index f8abf13b8..dba170687 100644 --- a/internal/ceres/cuda_buffer.h +++ b/internal/ceres/cuda_buffer.h @@ -66,8 +66,9 @@ class CudaBuffer { if (data_ != nullptr) { CHECK_EQ(cudaFree(data_), cudaSuccess); } - CHECK_EQ(cudaMalloc(&data_, size * sizeof(T)), cudaSuccess) << - "Failed to allocate " << size * sizeof(T) << " bytes of GPU memory"; + CHECK_EQ(cudaMalloc(&data_, size * sizeof(T)), cudaSuccess) + << "Failed to allocate " << size * sizeof(T) + << " bytes of GPU memory"; size_ = size; } } diff --git a/internal/ceres/cuda_kernels_test.cc b/internal/ceres/cuda_kernels_test.cc index fbdc4b389..83b922d12 100644 --- a/internal/ceres/cuda_kernels_test.cc +++ b/internal/ceres/cuda_kernels_test.cc @@ -52,10 +52,8 @@ TEST(CudaFP64ToFP32, SimpleConversions) { fp64_gpu.CopyFromCpuVector(fp64_cpu, cudaStreamDefault); CudaBuffer fp32_gpu; fp32_gpu.Reserve(fp64_cpu.size()); - CudaFP64ToFP32(fp64_gpu.data(), - fp32_gpu.data(), - fp64_cpu.size(), - cudaStreamDefault); + CudaFP64ToFP32( + fp64_gpu.data(), fp32_gpu.data(), fp64_cpu.size(), cudaStreamDefault); std::vector fp32_cpu(fp64_cpu.size()); fp32_gpu.CopyToCpu(fp32_cpu.data(), fp32_cpu.size()); for (int i = 0; i < fp32_cpu.size(); ++i) { @@ -65,11 +63,7 @@ TEST(CudaFP64ToFP32, SimpleConversions) { TEST(CudaFP64ToFP32, NumericallyExtremeValues) { std::vector fp64_cpu = { - DBL_MIN, - 10.0 * DBL_MIN, - DBL_MAX, - 0.1 * DBL_MAX - }; + DBL_MIN, 10.0 * DBL_MIN, DBL_MAX, 0.1 * DBL_MAX}; // First just make sure that the compiler has represented these values // accurately as fp64. EXPECT_GT(fp64_cpu[0], 0.0); @@ -80,10 +74,8 @@ TEST(CudaFP64ToFP32, NumericallyExtremeValues) { fp64_gpu.CopyFromCpuVector(fp64_cpu, cudaStreamDefault); CudaBuffer fp32_gpu; fp32_gpu.Reserve(fp64_cpu.size()); - CudaFP64ToFP32(fp64_gpu.data(), - fp32_gpu.data(), - fp64_cpu.size(), - cudaStreamDefault); + CudaFP64ToFP32( + fp64_gpu.data(), fp32_gpu.data(), fp64_cpu.size(), cudaStreamDefault); std::vector fp32_cpu(fp64_cpu.size()); fp32_gpu.CopyToCpu(fp32_cpu.data(), fp32_cpu.size()); EXPECT_EQ(fp32_cpu[0], 0.0f); @@ -98,10 +90,8 @@ TEST(CudaFP32ToFP64, SimpleConversions) { fp32_gpu.CopyFromCpuVector(fp32_cpu, cudaStreamDefault); CudaBuffer fp64_gpu; fp64_gpu.Reserve(fp32_cpu.size()); - CudaFP32ToFP64(fp32_gpu.data(), - fp64_gpu.data(), - fp32_cpu.size(), - cudaStreamDefault); + CudaFP32ToFP64( + fp32_gpu.data(), fp64_gpu.data(), fp32_cpu.size(), cudaStreamDefault); std::vector fp64_cpu(fp32_cpu.size()); fp64_gpu.CopyToCpu(fp64_cpu.data(), fp64_cpu.size()); for (int i = 0; i < fp64_cpu.size(); ++i) { @@ -135,8 +125,8 @@ TEST(CudaSetZeroFP64, NonZeroInput) { TEST(CudaDsxpy, DoubleValues) { std::vector fp32_cpu_a = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}; - std::vector fp64_cpu_b = - {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}; + std::vector fp64_cpu_b = { + 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}; CudaBuffer fp32_gpu_a; fp32_gpu_a.CopyFromCpuVector(fp32_cpu_a, cudaStreamDefault); CudaBuffer fp64_gpu_b; diff --git a/internal/ceres/dense_sparse_matrix.cc b/internal/ceres/dense_sparse_matrix.cc index e71546b86..67c3d2b59 100644 --- a/internal/ceres/dense_sparse_matrix.cc +++ b/internal/ceres/dense_sparse_matrix.cc @@ -59,11 +59,13 @@ DenseSparseMatrix::DenseSparseMatrix(Matrix m) : m_(std::move(m)) {} void DenseSparseMatrix::SetZero() { m_.setZero(); } -void DenseSparseMatrix::RightMultiply(const double* x, double* y) const { +void DenseSparseMatrix::RightMultiplyAndAccumulate(const double* x, + double* y) const { VectorRef(y, num_rows()).noalias() += m_ * ConstVectorRef(x, num_cols()); } -void DenseSparseMatrix::LeftMultiply(const double* x, double* y) const { +void DenseSparseMatrix::LeftMultiplyAndAccumulate(const double* x, + double* y) const { VectorRef(y, num_cols()).noalias() += m_.transpose() * ConstVectorRef(x, num_rows()); } diff --git a/internal/ceres/dense_sparse_matrix.h b/internal/ceres/dense_sparse_matrix.h index 160a59133..5fc71c20d 100644 --- a/internal/ceres/dense_sparse_matrix.h +++ b/internal/ceres/dense_sparse_matrix.h @@ -53,8 +53,8 @@ class CERES_NO_EXPORT DenseSparseMatrix final : public SparseMatrix { // SparseMatrix interface. void SetZero() final; - void RightMultiply(const double* x, double* y) const final; - void LeftMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; + void LeftMultiplyAndAccumulate(const double* x, double* y) const final; void SquaredColumnNorm(double* x) const final; void ScaleColumns(const double* scale) final; void ToDenseMatrix(Matrix* dense_matrix) const final; diff --git a/internal/ceres/dense_sparse_matrix_test.cc b/internal/ceres/dense_sparse_matrix_test.cc index 6bce0b437..0ba9fc235 100644 --- a/internal/ceres/dense_sparse_matrix_test.cc +++ b/internal/ceres/dense_sparse_matrix_test.cc @@ -59,8 +59,8 @@ static void CompareMatrices(const SparseMatrix* a, const SparseMatrix* b) { Vector y_a = Vector::Zero(num_rows); Vector y_b = Vector::Zero(num_rows); - a->RightMultiply(x.data(), y_a.data()); - b->RightMultiply(x.data(), y_b.data()); + a->RightMultiplyAndAccumulate(x.data(), y_a.data()); + b->RightMultiplyAndAccumulate(x.data(), y_b.data()); EXPECT_EQ((y_a - y_b).norm(), 0); } @@ -88,7 +88,7 @@ class DenseSparseMatrixTest : public ::testing::Test { std::unique_ptr dsm; }; -TEST_F(DenseSparseMatrixTest, RightMultiply) { +TEST_F(DenseSparseMatrixTest, RightMultiplyAndAccumulate) { CompareMatrices(tsm.get(), dsm.get()); // Try with a not entirely zero vector to verify column interactions, which @@ -100,13 +100,13 @@ TEST_F(DenseSparseMatrixTest, RightMultiply) { Vector b1 = Vector::Zero(num_rows); Vector b2 = Vector::Zero(num_rows); - tsm->RightMultiply(a.data(), b1.data()); - dsm->RightMultiply(a.data(), b2.data()); + tsm->RightMultiplyAndAccumulate(a.data(), b1.data()); + dsm->RightMultiplyAndAccumulate(a.data(), b2.data()); EXPECT_EQ((b1 - b2).norm(), 0); } -TEST_F(DenseSparseMatrixTest, LeftMultiply) { +TEST_F(DenseSparseMatrixTest, LeftMultiplyAndAccumulate) { for (int i = 0; i < num_rows; ++i) { Vector a = Vector::Zero(num_rows); a(i) = 1.0; @@ -114,8 +114,8 @@ TEST_F(DenseSparseMatrixTest, LeftMultiply) { Vector b1 = Vector::Zero(num_cols); Vector b2 = Vector::Zero(num_cols); - tsm->LeftMultiply(a.data(), b1.data()); - dsm->LeftMultiply(a.data(), b2.data()); + tsm->LeftMultiplyAndAccumulate(a.data(), b1.data()); + dsm->LeftMultiplyAndAccumulate(a.data(), b2.data()); EXPECT_EQ((b1 - b2).norm(), 0); } @@ -129,8 +129,8 @@ TEST_F(DenseSparseMatrixTest, LeftMultiply) { Vector b1 = Vector::Zero(num_cols); Vector b2 = Vector::Zero(num_cols); - tsm->LeftMultiply(a.data(), b1.data()); - dsm->LeftMultiply(a.data(), b2.data()); + tsm->LeftMultiplyAndAccumulate(a.data(), b1.data()); + dsm->LeftMultiplyAndAccumulate(a.data(), b2.data()); EXPECT_EQ((b1 - b2).norm(), 0); } diff --git a/internal/ceres/dogleg_strategy.cc b/internal/ceres/dogleg_strategy.cc index ac8c7d779..0db57de1c 100644 --- a/internal/ceres/dogleg_strategy.cc +++ b/internal/ceres/dogleg_strategy.cc @@ -175,7 +175,7 @@ TrustRegionStrategy::Summary DoglegStrategy::ComputeStep( void DoglegStrategy::ComputeGradient(SparseMatrix* jacobian, const double* residuals) { gradient_.setZero(); - jacobian->LeftMultiply(residuals, gradient_.data()); + jacobian->LeftMultiplyAndAccumulate(residuals, gradient_.data()); gradient_.array() /= diagonal_.array(); } @@ -188,7 +188,7 @@ void DoglegStrategy::ComputeCauchyPoint(SparseMatrix* jacobian) { // The Jacobian is scaled implicitly by computing J * (D^-1 * (D^-1 * g)) // instead of (J * D^-1) * (D^-1 * g). Vector scaled_gradient = (gradient_.array() / diagonal_.array()).matrix(); - jacobian->RightMultiply(scaled_gradient.data(), Jg.data()); + jacobian->RightMultiplyAndAccumulate(scaled_gradient.data(), Jg.data()); alpha_ = gradient_.squaredNorm() / Jg.squaredNorm(); } @@ -706,9 +706,9 @@ bool DoglegStrategy::ComputeSubspaceModel(SparseMatrix* jacobian) { Vector tmp; tmp = (subspace_basis_.col(0).array() / diagonal_.array()).matrix(); - jacobian->RightMultiply(tmp.data(), Jb.row(0).data()); + jacobian->RightMultiplyAndAccumulate(tmp.data(), Jb.row(0).data()); tmp = (subspace_basis_.col(1).array() / diagonal_.array()).matrix(); - jacobian->RightMultiply(tmp.data(), Jb.row(1).data()); + jacobian->RightMultiplyAndAccumulate(tmp.data(), Jb.row(1).data()); subspace_B_ = Jb * Jb.transpose(); diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc index 992d48c4a..81cf933a0 100644 --- a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc +++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc @@ -64,7 +64,7 @@ LinearSolver::Summary DynamicSparseNormalCholeskySolver::SolveImpl( double* x) { const int num_cols = A->num_cols(); VectorRef(x, num_cols).setZero(); - A->LeftMultiply(b, x); + A->LeftMultiplyAndAccumulate(b, x); if (per_solve_options.D != nullptr) { // Temporarily append a diagonal block to the A matrix, but undo diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc index f9ff44353..8d66022ff 100644 --- a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc +++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc @@ -72,7 +72,7 @@ class DynamicSparseNormalCholeskySolverTest : public ::testing::Test { Vector rhs(A_->num_cols()); rhs.setZero(); - A_->LeftMultiply(b_.get(), rhs.data()); + A_->LeftMultiplyAndAccumulate(b_.get(), rhs.data()); Vector expected_solution = lhs.llt().solve(rhs); std::unique_ptr solver(LinearSolver::Create(options)); diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc index 7946a568b..751612f9c 100644 --- a/internal/ceres/implicit_schur_complement.cc +++ b/internal/ceres/implicit_schur_complement.cc @@ -96,23 +96,24 @@ void ImplicitSchurComplement::Init(const BlockSparseMatrix& A, // By breaking it down into individual matrix vector products // involving the matrices E and F. This is implemented using a // PartitionedMatrixView of the input matrix A. -void ImplicitSchurComplement::RightMultiply(const double* x, double* y) const { +void ImplicitSchurComplement::RightMultiplyAndAccumulate(const double* x, + double* y) const { // y1 = F x tmp_rows_.setZero(); - A_->RightMultiplyF(x, tmp_rows_.data()); + A_->RightMultiplyAndAccumulateF(x, tmp_rows_.data()); // y2 = E' y1 tmp_e_cols_.setZero(); - A_->LeftMultiplyE(tmp_rows_.data(), tmp_e_cols_.data()); + A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(), tmp_e_cols_.data()); // y3 = -(E'E)^-1 y2 tmp_e_cols_2_.setZero(); - block_diagonal_EtE_inverse_->RightMultiply(tmp_e_cols_.data(), - tmp_e_cols_2_.data()); + block_diagonal_EtE_inverse_->RightMultiplyAndAccumulate(tmp_e_cols_.data(), + tmp_e_cols_2_.data()); tmp_e_cols_2_ *= -1.0; // y1 = y1 + E y3 - A_->RightMultiplyE(tmp_e_cols_2_.data(), tmp_rows_.data()); + A_->RightMultiplyAndAccumulateE(tmp_e_cols_2_.data(), tmp_rows_.data()); // y5 = D * x if (D_ != nullptr) { @@ -125,7 +126,7 @@ void ImplicitSchurComplement::RightMultiply(const double* x, double* y) const { } // y = y5 + F' y1 - A_->LeftMultiplyF(tmp_rows_.data(), y); + A_->LeftMultiplyAndAccumulateF(tmp_rows_.data(), y); } // Given a block diagonal matrix and an optional array of diagonal @@ -153,8 +154,8 @@ void ImplicitSchurComplement::AddDiagonalAndInvert( } } -// Similar to RightMultiply, use the block structure of the matrix A -// to compute y = (E'E)^-1 (E'b - E'F x). +// Similar to RightMultiplyAndAccumulate, use the block structure of the matrix +// A to compute y = (E'E)^-1 (E'b - E'F x). void ImplicitSchurComplement::BackSubstitute(const double* x, double* y) { const int num_cols_e = A_->num_cols_e(); const int num_cols_f = A_->num_cols_f(); @@ -163,18 +164,19 @@ void ImplicitSchurComplement::BackSubstitute(const double* x, double* y) { // y1 = F x tmp_rows_.setZero(); - A_->RightMultiplyF(x, tmp_rows_.data()); + A_->RightMultiplyAndAccumulateF(x, tmp_rows_.data()); // y2 = b - y1 tmp_rows_ = ConstVectorRef(b_, num_rows) - tmp_rows_; // y3 = E' y2 tmp_e_cols_.setZero(); - A_->LeftMultiplyE(tmp_rows_.data(), tmp_e_cols_.data()); + A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(), tmp_e_cols_.data()); // y = (E'E)^-1 y3 VectorRef(y, num_cols).setZero(); - block_diagonal_EtE_inverse_->RightMultiply(tmp_e_cols_.data(), y); + block_diagonal_EtE_inverse_->RightMultiplyAndAccumulate(tmp_e_cols_.data(), + y); // The full solution vector y has two blocks. The first block of // variables corresponds to the eliminated variables, which we just @@ -193,22 +195,23 @@ void ImplicitSchurComplement::BackSubstitute(const double* x, double* y) { void ImplicitSchurComplement::UpdateRhs() { // y1 = E'b tmp_e_cols_.setZero(); - A_->LeftMultiplyE(b_, tmp_e_cols_.data()); + A_->LeftMultiplyAndAccumulateE(b_, tmp_e_cols_.data()); // y2 = (E'E)^-1 y1 Vector y2 = Vector::Zero(A_->num_cols_e()); - block_diagonal_EtE_inverse_->RightMultiply(tmp_e_cols_.data(), y2.data()); + block_diagonal_EtE_inverse_->RightMultiplyAndAccumulate(tmp_e_cols_.data(), + y2.data()); // y3 = E y2 tmp_rows_.setZero(); - A_->RightMultiplyE(y2.data(), tmp_rows_.data()); + A_->RightMultiplyAndAccumulateE(y2.data(), tmp_rows_.data()); // y3 = b - y3 tmp_rows_ = ConstVectorRef(b_, A_->num_rows()) - tmp_rows_; // rhs = F' y3 rhs_.setZero(); - A_->LeftMultiplyF(tmp_rows_.data(), rhs_.data()); + A_->LeftMultiplyAndAccumulateF(tmp_rows_.data(), rhs_.data()); } } // namespace ceres::internal diff --git a/internal/ceres/implicit_schur_complement.h b/internal/ceres/implicit_schur_complement.h index 75e05a42e..8fcc309ee 100644 --- a/internal/ceres/implicit_schur_complement.h +++ b/internal/ceres/implicit_schur_complement.h @@ -85,9 +85,9 @@ class BlockSparseMatrix; // complement using the PartitionedMatrixView object. // // THREAD SAFETY: This class is not thread safe. In particular, the -// RightMultiply (and the LeftMultiply) methods are not thread safe as -// they depend on mutable arrays used for the temporaries needed to -// compute the product y += Sx; +// RightMultiplyAndAccumulate (and the LeftMultiplyAndAccumulate) methods are +// not thread safe as they depend on mutable arrays used for the temporaries +// needed to compute the product y += Sx; class CERES_NO_EXPORT ImplicitSchurComplement final : public LinearOperator { public: // num_eliminate_blocks is the number of E blocks in the matrix @@ -114,12 +114,12 @@ class CERES_NO_EXPORT ImplicitSchurComplement final : public LinearOperator { void Init(const BlockSparseMatrix& A, const double* D, const double* b); // y += Sx, where S is the Schur complement. - void RightMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; // The Schur complement is a symmetric positive definite matrix, // thus the left and right multiply operators are the same. - void LeftMultiply(const double* x, double* y) const final { - RightMultiply(x, y); + void LeftMultiplyAndAccumulate(const double* x, double* y) const final { + RightMultiplyAndAccumulate(x, y); } // y = (E'E)^-1 (E'b - E'F x). Given an estimate of the solution to @@ -155,7 +155,7 @@ class CERES_NO_EXPORT ImplicitSchurComplement final : public LinearOperator { Vector rhs_; - // Temporary storage vectors used to implement RightMultiply. + // Temporary storage vectors used to implement RightMultiplyAndAccumulate. mutable Vector tmp_rows_; mutable Vector tmp_e_cols_; mutable Vector tmp_e_cols_2_; diff --git a/internal/ceres/implicit_schur_complement_test.cc b/internal/ceres/implicit_schur_complement_test.cc index baa381afb..0ebde31e9 100644 --- a/internal/ceres/implicit_schur_complement_test.cc +++ b/internal/ceres/implicit_schur_complement_test.cc @@ -147,7 +147,7 @@ class ImplicitSchurComplementTest : public ::testing::Test { y = lhs * x; Vector z(num_sc_cols); - isc.RightMultiply(x.data(), z.data()); + isc.RightMultiplyAndAccumulate(x.data(), z.data()); // The i^th column of the implicit schur complement is the same as // the explicit schur complement. diff --git a/internal/ceres/iterative_refiner.cc b/internal/ceres/iterative_refiner.cc index 90ff5114a..aaeefa332 100644 --- a/internal/ceres/iterative_refiner.cc +++ b/internal/ceres/iterative_refiner.cc @@ -62,7 +62,7 @@ void SparseIterativeRefiner::Refine(const SparseMatrix& lhs, for (int i = 0; i < max_num_iterations_; ++i) { // residual = rhs - lhs * solution lhs_x_solution_.setZero(); - lhs.RightMultiply(solution_ptr, lhs_x_solution_.data()); + lhs.RightMultiplyAndAccumulate(solution_ptr, lhs_x_solution_.data()); residual_ = rhs - lhs_x_solution_; // solution += lhs^-1 residual cholesky->Solve(residual_.data(), correction_.data(), &ignored_message); diff --git a/internal/ceres/iterative_refiner_test.cc b/internal/ceres/iterative_refiner_test.cc index 9ea2340e5..49e379e92 100644 --- a/internal/ceres/iterative_refiner_test.cc +++ b/internal/ceres/iterative_refiner_test.cc @@ -58,13 +58,13 @@ class FakeSparseMatrix : public SparseMatrix { explicit FakeSparseMatrix(Matrix m) : m_(std::move(m)) {} // y += Ax - void RightMultiply(const double* x, double* y) const final { + void RightMultiplyAndAccumulate(const double* x, double* y) const final { VectorRef(y, m_.cols()) += m_ * ConstVectorRef(x, m_.cols()); } // y += A'x - void LeftMultiply(const double* x, double* y) const final { + void LeftMultiplyAndAccumulate(const double* x, double* y) const final { // We will assume that this is a symmetric matrix. - RightMultiply(x, y); + RightMultiplyAndAccumulate(x, y); } double* mutable_values() final { return m_.data(); } diff --git a/internal/ceres/line_search_direction.cc b/internal/ceres/line_search_direction.cc index e93d3e907..f14292fba 100644 --- a/internal/ceres/line_search_direction.cc +++ b/internal/ceres/line_search_direction.cc @@ -120,8 +120,8 @@ class CERES_NO_EXPORT LBFGS final : public LineSearchDirection { current.gradient - previous.gradient); search_direction->setZero(); - low_rank_inverse_hessian_.RightMultiply(current.gradient.data(), - search_direction->data()); + low_rank_inverse_hessian_.RightMultiplyAndAccumulate( + current.gradient.data(), search_direction->data()); *search_direction *= -1.0; if (search_direction->dot(current.gradient) >= 0.0) { diff --git a/internal/ceres/linear_operator.h b/internal/ceres/linear_operator.h index cab87e7aa..8a1e902d2 100644 --- a/internal/ceres/linear_operator.h +++ b/internal/ceres/linear_operator.h @@ -46,16 +46,16 @@ class CERES_NO_EXPORT LinearOperator { virtual ~LinearOperator(); // y = y + Ax; - virtual void RightMultiply(const double* x, double* y) const = 0; + virtual void RightMultiplyAndAccumulate(const double* x, double* y) const = 0; // y = y + A'x; - virtual void LeftMultiply(const double* x, double* y) const = 0; + virtual void LeftMultiplyAndAccumulate(const double* x, double* y) const = 0; - virtual void RightMultiply(const Vector& x, Vector& y) const { - RightMultiply(x.data(), y.data()); + virtual void RightMultiplyAndAccumulate(const Vector& x, Vector& y) const { + RightMultiplyAndAccumulate(x.data(), y.data()); } - virtual void LeftMultiply(const Vector& x, Vector& y) const { - LeftMultiply(x.data(), y.data()); + virtual void LeftMultiplyAndAccumulate(const Vector& x, Vector& y) const { + LeftMultiplyAndAccumulate(x.data(), y.data()); } virtual int num_rows() const = 0; diff --git a/internal/ceres/low_rank_inverse_hessian.cc b/internal/ceres/low_rank_inverse_hessian.cc index 42827e2ac..c47184486 100644 --- a/internal/ceres/low_rank_inverse_hessian.cc +++ b/internal/ceres/low_rank_inverse_hessian.cc @@ -116,8 +116,8 @@ bool LowRankInverseHessian::Update(const Vector& delta_x, return true; } -void LowRankInverseHessian::RightMultiply(const double* x_ptr, - double* y_ptr) const { +void LowRankInverseHessian::RightMultiplyAndAccumulate(const double* x_ptr, + double* y_ptr) const { ConstVectorRef gradient(x_ptr, num_parameters_); VectorRef search_direction(y_ptr, num_parameters_); diff --git a/internal/ceres/low_rank_inverse_hessian.h b/internal/ceres/low_rank_inverse_hessian.h index de30f544b..878db81d4 100644 --- a/internal/ceres/low_rank_inverse_hessian.h +++ b/internal/ceres/low_rank_inverse_hessian.h @@ -64,7 +64,7 @@ class CERES_NO_EXPORT LowRankInverseHessian final : public LinearOperator { // num_parameters is the row/column size of the Hessian. // max_num_corrections is the rank of the Hessian approximation. // use_approximate_eigenvalue_scaling controls whether the initial - // inverse Hessian used during Right/LeftMultiply() is scaled by + // inverse Hessian used during Right/LeftMultiplyAndAccumulate() is scaled by // the approximate eigenvalue of the true inverse Hessian at the // current operating point. // The approximation uses: @@ -83,9 +83,9 @@ class CERES_NO_EXPORT LowRankInverseHessian final : public LinearOperator { bool Update(const Vector& delta_x, const Vector& delta_gradient); // LinearOperator interface - void RightMultiply(const double* x, double* y) const final; - void LeftMultiply(const double* x, double* y) const final { - RightMultiply(x, y); + void RightMultiplyAndAccumulate(const double* x, double* y) const final; + void LeftMultiplyAndAccumulate(const double* x, double* y) const final { + RightMultiplyAndAccumulate(x, y); } int num_rows() const final { return num_parameters_; } int num_cols() const final { return num_parameters_; } diff --git a/internal/ceres/partitioned_matrix_view.h b/internal/ceres/partitioned_matrix_view.h index 130572035..7fb1a095d 100644 --- a/internal/ceres/partitioned_matrix_view.h +++ b/internal/ceres/partitioned_matrix_view.h @@ -67,16 +67,18 @@ class CERES_NO_EXPORT PartitionedMatrixViewBase { virtual ~PartitionedMatrixViewBase(); // y += E'x - virtual void LeftMultiplyE(const double* x, double* y) const = 0; + virtual void LeftMultiplyAndAccumulateE(const double* x, double* y) const = 0; // y += F'x - virtual void LeftMultiplyF(const double* x, double* y) const = 0; + virtual void LeftMultiplyAndAccumulateF(const double* x, double* y) const = 0; // y += Ex - virtual void RightMultiplyE(const double* x, double* y) const = 0; + virtual void RightMultiplyAndAccumulateE(const double* x, + double* y) const = 0; // y += Fx - virtual void RightMultiplyF(const double* x, double* y) const = 0; + virtual void RightMultiplyAndAccumulateF(const double* x, + double* y) const = 0; // Create and return the block diagonal of the matrix E'E. virtual std::unique_ptr CreateBlockDiagonalEtE() const = 0; @@ -124,10 +126,10 @@ class CERES_NO_EXPORT PartitionedMatrixView final // num_col_blocks_a column blocks. PartitionedMatrixView(const BlockSparseMatrix& matrix, int num_col_blocks_e); - void LeftMultiplyE(const double* x, double* y) const final; - void LeftMultiplyF(const double* x, double* y) const final; - void RightMultiplyE(const double* x, double* y) const final; - void RightMultiplyF(const double* x, double* y) const final; + void LeftMultiplyAndAccumulateE(const double* x, double* y) const final; + void LeftMultiplyAndAccumulateF(const double* x, double* y) const final; + void RightMultiplyAndAccumulateE(const double* x, double* y) const final; + void RightMultiplyAndAccumulateF(const double* x, double* y) const final; std::unique_ptr CreateBlockDiagonalEtE() const final; std::unique_ptr CreateBlockDiagonalFtF() const final; void UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const final; diff --git a/internal/ceres/partitioned_matrix_view_impl.h b/internal/ceres/partitioned_matrix_view_impl.h index a33b86d86..215066062 100644 --- a/internal/ceres/partitioned_matrix_view_impl.h +++ b/internal/ceres/partitioned_matrix_view_impl.h @@ -87,7 +87,7 @@ PartitionedMatrixView:: template void PartitionedMatrixView:: - RightMultiplyE(const double* x, double* y) const { + RightMultiplyAndAccumulateE(const double* x, double* y) const { const CompressedRowBlockStructure* bs = matrix_.block_structure(); // Iterate over the first num_row_blocks_e_ row blocks, and multiply @@ -111,7 +111,7 @@ void PartitionedMatrixView:: template void PartitionedMatrixView:: - RightMultiplyF(const double* x, double* y) const { + RightMultiplyAndAccumulateF(const double* x, double* y) const { const CompressedRowBlockStructure* bs = matrix_.block_structure(); // Iterate over row blocks, and if the row block is in E, then @@ -157,7 +157,7 @@ void PartitionedMatrixView:: template void PartitionedMatrixView:: - LeftMultiplyE(const double* x, double* y) const { + LeftMultiplyAndAccumulateE(const double* x, double* y) const { const CompressedRowBlockStructure* bs = matrix_.block_structure(); // Iterate over the first num_row_blocks_e_ row blocks, and multiply @@ -181,7 +181,7 @@ void PartitionedMatrixView:: template void PartitionedMatrixView:: - LeftMultiplyF(const double* x, double* y) const { + LeftMultiplyAndAccumulateF(const double* x, double* y) const { const CompressedRowBlockStructure* bs = matrix_.block_structure(); // Iterate over row blocks, and if the row block is in E, then @@ -289,8 +289,8 @@ PartitionedMatrixView:: return block_diagonal; } -// Similar to the code in RightMultiplyE, except instead of the matrix -// vector multiply its an outer product. +// Similar to the code in RightMultiplyAndAccumulateE, except instead of the +// matrix vector multiply its an outer product. // // block_diagonal = block_diagonal(E'E) // @@ -322,8 +322,8 @@ void PartitionedMatrixView:: } } -// Similar to the code in RightMultiplyF, except instead of the matrix -// vector multiply its an outer product. +// Similar to the code in RightMultiplyAndAccumulateF, except instead of the +// matrix vector multiply its an outer product. // // block_diagonal = block_diagonal(F'F) // diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc index e43e32f78..cb3dd141f 100644 --- a/internal/ceres/partitioned_matrix_view_test.cc +++ b/internal/ceres/partitioned_matrix_view_test.cc @@ -85,7 +85,7 @@ TEST_F(PartitionedMatrixViewTest, DimensionsTest) { EXPECT_EQ(pmv_->num_rows(), A_->num_rows()); } -TEST_F(PartitionedMatrixViewTest, RightMultiplyE) { +TEST_F(PartitionedMatrixViewTest, RightMultiplyAndAccumulateE) { Vector x1(pmv_->num_cols_e()); Vector x2(pmv_->num_cols()); x2.setZero(); @@ -95,17 +95,17 @@ TEST_F(PartitionedMatrixViewTest, RightMultiplyE) { } Vector y1 = Vector::Zero(pmv_->num_rows()); - pmv_->RightMultiplyE(x1.data(), y1.data()); + pmv_->RightMultiplyAndAccumulateE(x1.data(), y1.data()); Vector y2 = Vector::Zero(pmv_->num_rows()); - A_->RightMultiply(x2.data(), y2.data()); + A_->RightMultiplyAndAccumulate(x2.data(), y2.data()); for (int i = 0; i < pmv_->num_rows(); ++i) { EXPECT_NEAR(y1(i), y2(i), kEpsilon); } } -TEST_F(PartitionedMatrixViewTest, RightMultiplyF) { +TEST_F(PartitionedMatrixViewTest, RightMultiplyAndAccumulateF) { Vector x1(pmv_->num_cols_f()); Vector x2 = Vector::Zero(pmv_->num_cols()); @@ -115,17 +115,17 @@ TEST_F(PartitionedMatrixViewTest, RightMultiplyF) { } Vector y1 = Vector::Zero(pmv_->num_rows()); - pmv_->RightMultiplyF(x1.data(), y1.data()); + pmv_->RightMultiplyAndAccumulateF(x1.data(), y1.data()); Vector y2 = Vector::Zero(pmv_->num_rows()); - A_->RightMultiply(x2.data(), y2.data()); + A_->RightMultiplyAndAccumulate(x2.data(), y2.data()); for (int i = 0; i < pmv_->num_rows(); ++i) { EXPECT_NEAR(y1(i), y2(i), kEpsilon); } } -TEST_F(PartitionedMatrixViewTest, LeftMultiply) { +TEST_F(PartitionedMatrixViewTest, LeftMultiplyAndAccumulate) { Vector x = Vector::Zero(pmv_->num_rows()); for (int i = 0; i < pmv_->num_rows(); ++i) { x(i) = RandDouble(); @@ -135,9 +135,9 @@ TEST_F(PartitionedMatrixViewTest, LeftMultiply) { Vector y1 = Vector::Zero(pmv_->num_cols_e()); Vector y2 = Vector::Zero(pmv_->num_cols_f()); - A_->LeftMultiply(x.data(), y.data()); - pmv_->LeftMultiplyE(x.data(), y1.data()); - pmv_->LeftMultiplyF(x.data(), y2.data()); + A_->LeftMultiplyAndAccumulate(x.data(), y.data()); + pmv_->LeftMultiplyAndAccumulateE(x.data(), y1.data()); + pmv_->LeftMultiplyAndAccumulateF(x.data(), y2.data()); for (int i = 0; i < pmv_->num_cols(); ++i) { EXPECT_NEAR(y(i), diff --git a/internal/ceres/preconditioner.cc b/internal/ceres/preconditioner.cc index 391e1b5bc..33174dea5 100644 --- a/internal/ceres/preconditioner.cc +++ b/internal/ceres/preconditioner.cc @@ -60,9 +60,9 @@ bool SparseMatrixPreconditionerWrapper::UpdateImpl(const SparseMatrix& A, return true; } -void SparseMatrixPreconditionerWrapper::RightMultiply(const double* x, - double* y) const { - matrix_->RightMultiply(x, y); +void SparseMatrixPreconditionerWrapper::RightMultiplyAndAccumulate( + const double* x, double* y) const { + matrix_->RightMultiplyAndAccumulate(x, y); } int SparseMatrixPreconditionerWrapper::num_rows() const { diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h index 75613fb50..2d343bd43 100644 --- a/internal/ceres/preconditioner.h +++ b/internal/ceres/preconditioner.h @@ -130,12 +130,13 @@ class CERES_NO_EXPORT Preconditioner : public LinearOperator { virtual bool Update(const LinearOperator& A, const double* D) = 0; // LinearOperator interface. Since the operator is symmetric, - // LeftMultiply and num_cols are just calls to RightMultiply and - // num_rows respectively. Update() must be called before - // RightMultiply can be called. - void RightMultiply(const double* x, double* y) const override = 0; - void LeftMultiply(const double* x, double* y) const override { - return RightMultiply(x, y); + // LeftMultiplyAndAccumulate and num_cols are just calls to + // RightMultiplyAndAccumulate and num_rows respectively. Update() must be + // called before RightMultiplyAndAccumulate can be called. + void RightMultiplyAndAccumulate(const double* x, + double* y) const override = 0; + void LeftMultiplyAndAccumulate(const double* x, double* y) const override { + return RightMultiplyAndAccumulate(x, y); } int num_rows() const override = 0; @@ -148,7 +149,7 @@ class CERES_NO_EXPORT IdentityPreconditioner : public Preconditioner { bool Update(const LinearOperator& A, const double* D) final { return true; } - void RightMultiply(const double* x, double* y) const final { + void RightMultiplyAndAccumulate(const double* x, double* y) const final { VectorRef(y, num_rows_) += ConstVectorRef(x, num_rows_); } @@ -189,7 +190,7 @@ class CERES_NO_EXPORT SparseMatrixPreconditionerWrapper final ~SparseMatrixPreconditionerWrapper() override; // Preconditioner interface - void RightMultiply(const double* x, double* y) const override; + void RightMultiplyAndAccumulate(const double* x, double* y) const override; int num_rows() const override; private: diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index 4e478373e..1f10ac27a 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc @@ -70,8 +70,8 @@ class BlockRandomAccessSparseMatrixAdapter virtual ~BlockRandomAccessSparseMatrixAdapter() final {} - void RightMultiply(const Vector& x, Vector& y) final { - m_.SymmetricRightMultiply(x.data(), y.data()); + void RightMultiplyAndAccumulate(const Vector& x, Vector& y) final { + m_.SymmetricRightMultiplyAndAccumulate(x.data(), y.data()); } private: @@ -88,8 +88,8 @@ class BlockRandomAccessDiagonalMatrixAdapter final virtual ~BlockRandomAccessDiagonalMatrixAdapter() final {} // y = y + Ax; - void RightMultiply(const Vector& x, Vector& y) final { - m_.RightMultiply(x.data(), y.data()); + void RightMultiplyAndAccumulate(const Vector& x, Vector& y) final { + m_.RightMultiplyAndAccumulate(x.data(), y.data()); } private: diff --git a/internal/ceres/schur_jacobi_preconditioner.cc b/internal/ceres/schur_jacobi_preconditioner.cc index d452ba496..331792763 100644 --- a/internal/ceres/schur_jacobi_preconditioner.cc +++ b/internal/ceres/schur_jacobi_preconditioner.cc @@ -91,9 +91,9 @@ bool SchurJacobiPreconditioner::UpdateImpl(const BlockSparseMatrix& A, return true; } -void SchurJacobiPreconditioner::RightMultiply(const double* x, - double* y) const { - m_->RightMultiply(x, y); +void SchurJacobiPreconditioner::RightMultiplyAndAccumulate(const double* x, + double* y) const { + m_->RightMultiplyAndAccumulate(x, y); } int SchurJacobiPreconditioner::num_rows() const { return m_->num_rows(); } diff --git a/internal/ceres/schur_jacobi_preconditioner.h b/internal/ceres/schur_jacobi_preconditioner.h index 76a7b3d8d..ddf471c13 100644 --- a/internal/ceres/schur_jacobi_preconditioner.h +++ b/internal/ceres/schur_jacobi_preconditioner.h @@ -71,7 +71,7 @@ class SchurEliminatorBase; // SchurJacobiPreconditioner preconditioner( // *A.block_structure(), options); // preconditioner.Update(A, nullptr); -// preconditioner.RightMultiply(x, y); +// preconditioner.RightMultiplyAndAccumulate(x, y); // class CERES_NO_EXPORT SchurJacobiPreconditioner : public BlockSparseMatrixPreconditioner { @@ -90,7 +90,7 @@ class CERES_NO_EXPORT SchurJacobiPreconditioner ~SchurJacobiPreconditioner() override; // Preconditioner interface. - void RightMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; int num_rows() const final; private: diff --git a/internal/ceres/sparse_matrix.h b/internal/ceres/sparse_matrix.h index 9fe33a30f..da6af184e 100644 --- a/internal/ceres/sparse_matrix.h +++ b/internal/ceres/sparse_matrix.h @@ -68,9 +68,10 @@ class CERES_NO_EXPORT SparseMatrix : public LinearOperator { ~SparseMatrix() override; // y += Ax; - void RightMultiply(const double* x, double* y) const override = 0; + void RightMultiplyAndAccumulate(const double* x, + double* y) const override = 0; // y += A'x; - void LeftMultiply(const double* x, double* y) const override = 0; + void LeftMultiplyAndAccumulate(const double* x, double* y) const override = 0; // In MATLAB notation sum(A.*A, 1) virtual void SquaredColumnNorm(double* x) const = 0; diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc index 949991b89..99205fa2c 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/internal/ceres/sparse_normal_cholesky_solver.cc @@ -71,7 +71,7 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl( xref.setZero(); rhs_.resize(num_cols); rhs_.setZero(); - A->LeftMultiply(b, rhs_.data()); + A->LeftMultiplyAndAccumulate(b, rhs_.data()); event_logger.AddEvent("Compute RHS"); if (per_solve_options.D != nullptr) { diff --git a/internal/ceres/sparse_normal_cholesky_solver_test.cc b/internal/ceres/sparse_normal_cholesky_solver_test.cc index b7d4a3965..002b70761 100644 --- a/internal/ceres/sparse_normal_cholesky_solver_test.cc +++ b/internal/ceres/sparse_normal_cholesky_solver_test.cc @@ -75,7 +75,7 @@ class SparseNormalCholeskySolverTest : public ::testing::Test { Vector rhs(A_->num_cols()); rhs.setZero(); - A_->LeftMultiply(b_.get(), rhs.data()); + A_->LeftMultiplyAndAccumulate(b_.get(), rhs.data()); Vector expected_solution = lhs.llt().solve(rhs); std::unique_ptr solver(LinearSolver::Create(options)); diff --git a/internal/ceres/subset_preconditioner.cc b/internal/ceres/subset_preconditioner.cc index a22854599..5dc364f52 100644 --- a/internal/ceres/subset_preconditioner.cc +++ b/internal/ceres/subset_preconditioner.cc @@ -57,7 +57,8 @@ SubsetPreconditioner::SubsetPreconditioner(Preconditioner::Options options, SubsetPreconditioner::~SubsetPreconditioner() = default; -void SubsetPreconditioner::RightMultiply(const double* x, double* y) const { +void SubsetPreconditioner::RightMultiplyAndAccumulate(const double* x, + double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); std::string message; diff --git a/internal/ceres/subset_preconditioner.h b/internal/ceres/subset_preconditioner.h index 1f1c6ec3a..7139ca650 100644 --- a/internal/ceres/subset_preconditioner.h +++ b/internal/ceres/subset_preconditioner.h @@ -75,7 +75,7 @@ class CERES_NO_EXPORT SubsetPreconditioner ~SubsetPreconditioner() override; // Preconditioner interface - void RightMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; int num_rows() const final { return num_cols_; } int num_cols() const final { return num_cols_; } diff --git a/internal/ceres/subset_preconditioner_test.cc b/internal/ceres/subset_preconditioner_test.cc index 2f5e044f4..cd8669562 100644 --- a/internal/ceres/subset_preconditioner_test.cc +++ b/internal/ceres/subset_preconditioner_test.cc @@ -159,7 +159,7 @@ TEST_P(SubsetPreconditionerTest, foo) { EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected)); Vector actual(lhs->num_rows()); - preconditioner_->RightMultiply(rhs.data(), actual.data()); + preconditioner_->RightMultiplyAndAccumulate(rhs.data(), actual.data()); Matrix eigen_lhs; lhs->ToDenseMatrix(&eigen_lhs); diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc index b8f43adbd..49d367a1f 100644 --- a/internal/ceres/triplet_sparse_matrix.cc +++ b/internal/ceres/triplet_sparse_matrix.cc @@ -169,13 +169,15 @@ void TripletSparseMatrix::CopyData(const TripletSparseMatrix& orig) { } } -void TripletSparseMatrix::RightMultiply(const double* x, double* y) const { +void TripletSparseMatrix::RightMultiplyAndAccumulate(const double* x, + double* y) const { for (int i = 0; i < num_nonzeros_; ++i) { y[rows_[i]] += values_[i] * x[cols_[i]]; } } -void TripletSparseMatrix::LeftMultiply(const double* x, double* y) const { +void TripletSparseMatrix::LeftMultiplyAndAccumulate(const double* x, + double* y) const { for (int i = 0; i < num_nonzeros_; ++i) { y[cols_[i]] += values_[i] * x[rows_[i]]; } diff --git a/internal/ceres/triplet_sparse_matrix.h b/internal/ceres/triplet_sparse_matrix.h index a2ba1f198..c9624556d 100644 --- a/internal/ceres/triplet_sparse_matrix.h +++ b/internal/ceres/triplet_sparse_matrix.h @@ -65,8 +65,8 @@ class CERES_NO_EXPORT TripletSparseMatrix final : public SparseMatrix { // Implementation of the SparseMatrix interface. void SetZero() final; - void RightMultiply(const double* x, double* y) const final; - void LeftMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; + void LeftMultiplyAndAccumulate(const double* x, double* y) const final; void SquaredColumnNorm(double* x) const final; void ScaleColumns(const double* scale) final; void ToCRSMatrix(CRSMatrix* matrix) const; @@ -139,6 +139,7 @@ class CERES_NO_EXPORT TripletSparseMatrix final : public SparseMatrix { // 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 577d95990..b41d99180 100644 --- a/internal/ceres/triplet_sparse_matrix_test.cc +++ b/internal/ceres/triplet_sparse_matrix_test.cc @@ -28,11 +28,11 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) -#include "ceres/crs_matrix.h" #include "ceres/triplet_sparse_matrix.h" #include +#include "ceres/crs_matrix.h" #include "gtest/gtest.h" namespace ceres::internal { diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc index 739304ae7..6693e6eae 100644 --- a/internal/ceres/trust_region_minimizer.cc +++ b/internal/ceres/trust_region_minimizer.cc @@ -421,7 +421,8 @@ bool TrustRegionMinimizer::ComputeTrustRegionStep() { // = -f'J * step - step' * J' * J * step / 2 // = -(J * step)'(f + J * step / 2) model_residuals_.setZero(); - jacobian_->RightMultiply(trust_region_step_.data(), model_residuals_.data()); + jacobian_->RightMultiplyAndAccumulate(trust_region_step_.data(), + model_residuals_.data()); model_cost_change_ = -model_residuals_.dot(residuals_ + model_residuals_ / 2.0); diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc index 0a7b19f97..bdbe3c40a 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc @@ -422,8 +422,8 @@ LinearSolverTerminationType VisibilityBasedPreconditioner::Factorize() { return sparse_cholesky_->Factorize(lhs.get(), &message); } -void VisibilityBasedPreconditioner::RightMultiply(const double* x, - double* y) const { +void VisibilityBasedPreconditioner::RightMultiplyAndAccumulate( + const double* x, double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); CHECK(sparse_cholesky_ != nullptr); diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h index 39bdc9fea..f27eb9e5f 100644 --- a/internal/ceres/visibility_based_preconditioner.h +++ b/internal/ceres/visibility_based_preconditioner.h @@ -122,7 +122,7 @@ class SchurEliminatorBase; // VisibilityBasedPreconditioner preconditioner( // *A.block_structure(), options); // preconditioner.Update(A, nullptr); -// preconditioner.RightMultiply(x, y); +// preconditioner.RightMultiplyAndAccumulate(x, y); class CERES_NO_EXPORT VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner { public: @@ -140,7 +140,7 @@ class CERES_NO_EXPORT VisibilityBasedPreconditioner ~VisibilityBasedPreconditioner() override; // Preconditioner interface - void RightMultiply(const double* x, double* y) const final; + void RightMultiplyAndAccumulate(const double* x, double* y) const final; int num_rows() const final; friend class VisibilityBasedPreconditionerTest; diff --git a/internal/ceres/visibility_based_preconditioner_test.cc b/internal/ceres/visibility_based_preconditioner_test.cc index 8e0d5fd70..4cf1dbe7c 100644 --- a/internal/ceres/visibility_based_preconditioner_test.cc +++ b/internal/ceres/visibility_based_preconditioner_test.cc @@ -276,7 +276,7 @@ namespace ceres::internal { // y.setZero(); // z.setZero(); // x[i] = 1.0; -// preconditioner_->RightMultiply(x.data(), y.data()); +// preconditioner_->RightMultiplyAndAccumulate(x.data(), y.data()); // z = full_schur_complement // .selfadjointView() // .llt().solve(x);