Use solveInPlace correctly.

When solving a linear system using Eigen's dense Cholesky factorization if the right hand side of the linear system is the same vector that will store the solution, call solveInPlace instead of solve. Change-Id: I3e6d2f21ff420c25217cd87ee5d269fdfabbf19a
G-P-S · Apr 6, 2015 · 59d7f98 · 59d7f98
1 parent 10cbe85
commit 59d7f98
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 5 deletions.
diff --git a/internal/ceres/eigen_dense_cholesky.cc b/internal/ceres/eigen_dense_cholesky.cc
@@ -35,7 +35,7 @@ Eigen::ComputationInfo
 SolveUpperTriangularUsingCholesky(int size,
                                   const double* lhs_values,
                                   const double* rhs_values,
-                                  double* solution) {
+                                  double* solution_values) {
   ConstMatrixRef lhs(lhs_values, size, size);
 
   // On ARM we have experienced significant numerical problems with
@@ -51,8 +51,13 @@ SolveUpperTriangularUsingCholesky(int size,
 #endif
 
   if (cholesky.info() == Eigen::Success) {
-    ConstVectorRef rhs(rhs_values, size);
-    VectorRef(solution, size) = cholesky.solve(rhs);
+    VectorRef solution(solution_values, size);
+    if (solution_values == rhs_values) {
+      cholesky.solveInPlace(solution);
+    } else {
+      ConstVectorRef rhs(rhs_values, size);
+      solution = cholesky.solve(rhs);
+    }
   }
 
   return cholesky.info();

diff --git a/internal/ceres/eigen_dense_cholesky.h b/internal/ceres/eigen_dense_cholesky.h
@@ -57,8 +57,8 @@ InvertUpperTriangularUsingCholesky(int size,
 // factorization. rhs_values and solution can point to the same array.
 Eigen::ComputationInfo
 SolveUpperTriangularUsingCholesky(int size,
-                                  const double* lhs_values,
-                                  const double* rhs_values,
+                                  const double* lhs,
+                                  const double* rhs,
                                   double* solution);
 
 }  // namespace internal