1 files changed, 10 insertions, 76 deletions

diff --git a/extern/ceres/internal/ceres/dense_normal_cholesky_solver.cc b/extern/ceres/internal/ceres/dense_normal_cholesky_solver.cc
index 51c639097b6..30a0c023f51 100644
--- a/extern/ceres/internal/ceres/dense_normal_cholesky_solver.cc
+++ b/extern/ceres/internal/ceres/dense_normal_cholesky_solver.cc
@@ -30,13 +30,11 @@
 
 #include "ceres/dense_normal_cholesky_solver.h"
 
-#include <cstddef>
+#include <utility>
 
 #include "Eigen/Dense"
-#include "ceres/blas.h"
 #include "ceres/dense_sparse_matrix.h"
 #include "ceres/internal/eigen.h"
-#include "ceres/lapack.h"
 #include "ceres/linear_solver.h"
 #include "ceres/types.h"
 #include "ceres/wall_time.h"
@@ -45,32 +43,20 @@ namespace ceres {
 namespace internal {
 
 DenseNormalCholeskySolver::DenseNormalCholeskySolver(
-    const LinearSolver::Options& options)
-    : options_(options) {}
+    LinearSolver::Options options)
+    : options_(std::move(options)),
+      cholesky_(DenseCholesky::Create(options_)) {}
 
 LinearSolver::Summary DenseNormalCholeskySolver::SolveImpl(
     DenseSparseMatrix* A,
     const double* b,
     const LinearSolver::PerSolveOptions& per_solve_options,
     double* x) {
-  if (options_.dense_linear_algebra_library_type == EIGEN) {
-    return SolveUsingEigen(A, b, per_solve_options, x);
-  } else {
-    return SolveUsingLAPACK(A, b, per_solve_options, x);
-  }
-}
-
-LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen(
-    DenseSparseMatrix* A,
-    const double* b,
-    const LinearSolver::PerSolveOptions& per_solve_options,
-    double* x) {
   EventLogger event_logger("DenseNormalCholeskySolver::Solve");
 
   const int num_rows = A->num_rows();
   const int num_cols = A->num_cols();
 
-  ConstColMajorMatrixRef Aref = A->matrix();
   Matrix lhs(num_cols, num_cols);
   lhs.setZero();
 
@@ -81,12 +67,12 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen(
   // Using rankUpdate instead of GEMM, exposes the fact that its the
   // same matrix being multiplied with itself and that the product is
   // symmetric.
-  lhs.selfadjointView<Eigen::Upper>().rankUpdate(Aref.transpose());
+  lhs.selfadjointView<Eigen::Upper>().rankUpdate(A->matrix().transpose());
 
   //   rhs = A'b
-  Vector rhs = Aref.transpose() * ConstVectorRef(b, num_rows);
+  Vector rhs = A->matrix().transpose() * ConstVectorRef(b, num_rows);
 
-  if (per_solve_options.D != NULL) {
+  if (per_solve_options.D != nullptr) {
     ConstVectorRef D(per_solve_options.D, num_cols);
     lhs += D.array().square().matrix().asDiagonal();
   }
@@ -94,64 +80,12 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen(
 
   LinearSolver::Summary summary;
   summary.num_iterations = 1;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  Eigen::LLT<Matrix, Eigen::Upper> llt =
-      lhs.selfadjointView<Eigen::Upper>().llt();
-
-  if (llt.info() != Eigen::Success) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "Eigen LLT decomposition failed.";
-  } else {
-    summary.termination_type = LINEAR_SOLVER_SUCCESS;
-    summary.message = "Success.";
-  }
+  summary.termination_type = cholesky_->FactorAndSolve(
+      num_cols, lhs.data(), rhs.data(), x, &summary.message);
+  event_logger.AddEvent("FactorAndSolve");
 
-  VectorRef(x, num_cols) = llt.solve(rhs);
-  event_logger.AddEvent("Solve");
   return summary;
 }
 
-LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingLAPACK(
-    DenseSparseMatrix* A,
-    const double* b,
-    const LinearSolver::PerSolveOptions& per_solve_options,
-    double* x) {
-  EventLogger event_logger("DenseNormalCholeskySolver::Solve");
-
-  if (per_solve_options.D != NULL) {
-    // Temporarily append a diagonal block to the A matrix, but undo
-    // it before returning the matrix to the user.
-    A->AppendDiagonal(per_solve_options.D);
-  }
-
-  const int num_cols = A->num_cols();
-  Matrix lhs(num_cols, num_cols);
-  event_logger.AddEvent("Setup");
-
-  // lhs = A'A
-  //
-  // Note: This is a bit delicate, it assumes that the stride on this
-  // matrix is the same as the number of rows.
-  BLAS::SymmetricRankKUpdate(
-      A->num_rows(), num_cols, A->values(), true, 1.0, 0.0, lhs.data());
-
-  if (per_solve_options.D != NULL) {
-    // Undo the modifications to the matrix A.
-    A->RemoveDiagonal();
-  }
-
-  // TODO(sameeragarwal): Replace this with a gemv call for true blasness.
-  //   rhs = A'b
-  VectorRef(x, num_cols) =
-      A->matrix().transpose() * ConstVectorRef(b, A->num_rows());
-  event_logger.AddEvent("Product");
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 1;
-  summary.termination_type = LAPACK::SolveInPlaceUsingCholesky(
-      num_cols, lhs.data(), x, &summary.message);
-  event_logger.AddEvent("Solve");
-  return summary;
-}
 }  // namespace internal
 }  // namespace ceres