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<Eigen::Upper>() * 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<double>::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<Param> info) {
   return "UNSYMMETRIC";
 }
 
-class RightMultiplyTest : public ::testing::TestWithParam<Param> {};
+class RightMultiplyAndAccumulateTest : public ::testing::TestWithParam<Param> {
+};
 
-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<Param> {};
+class LeftMultiplyAndAccumulateTest : public ::testing::TestWithParam<Param> {};
 
-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 <typename DenseVectorType>
 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<Vector> {
   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<float> 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<float> 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<double> 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<float> 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<float> 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<double> 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<double> 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<float> fp32_cpu_a = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0};
-  std::vector<double> fp64_cpu_b =
-      {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0};
+  std::vector<double> fp64_cpu_b = {
+      1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0};
   CudaBuffer<float> fp32_gpu_a;
   fp32_gpu_a.CopyFromCpuVector(fp32_cpu_a, cudaStreamDefault);
   CudaBuffer<double> 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<DenseSparseMatrix> 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<LinearSolver> 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<BlockSparseMatrix> 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<BlockSparseMatrix> CreateBlockDiagonalEtE() const final;
   std::unique_ptr<BlockSparseMatrix> 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<kRowBlockSize, kEBlockSize, kFBlockSize>::
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    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<kRowBlockSize, kEBlockSize, kFBlockSize>::
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    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<kRowBlockSize, kEBlockSize, kFBlockSize>::
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    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<kRowBlockSize, kEBlockSize, kFBlockSize>::
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    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<kRowBlockSize, kEBlockSize, kFBlockSize>::
   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<kRowBlockSize, kEBlockSize, kFBlockSize>::
   }
 }
 
-// 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<LinearSolver> 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<TripletSparseMatrix> 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 <memory>
 
+#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<Eigen::Upper>()
 //         .llt().solve(x);