commit b43e73a03485f0fd0fe514e356ad8925731d3a81
author Sameer Agarwal <sameeragarwal@google.com> Fri Aug 01 12:09:09 2014 -0700
committer Sameer Agarwal <sameeragarwal@google.com> Fri Aug 01 12:10:04 2014 -0700
tree 8b570be7d441d859886fb0773075ad0b6ca233bb
parent 031598295c6b2f061c171b9b2338919f41b7eb0b

Simplify the Eigen code in SparseNormalCholeskySolver.

Simplifying some of the template handling, and remove the use
of SelfAdjointView as it is not needed. The solver itself takes
an argument for where the data is actually stored.

The performance of SparseNormalCholesky with EIGEN_SPARSE
seems to be on par with CX_SPARSE.

Change-Id: I69e22a144b447c052b6cbe59ef1aa33eae2dd9e3

internal/ceres/sparse_normal_cholesky_solver.cc

diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 4277480..3fd6d56 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -170,12 +170,6 @@
   // outer_product_ is a compressed row sparse matrix and in lower
   // triangular form, when mapped to a compressed column sparse
   // matrix, it becomes an upper triangular matrix.
-  //
-  // TODO(sameeragarwal): It is not clear to me if an upper triangular
-  // column major matrix is the way to go here, or if a lower
-  // triangular matrix is better. This will require some testing. If
-  // it turns out that the lower triangular is better, then the logic
-  // used to compute the outer product needs to be updated.
   Eigen::MappedSparseMatrix<double, Eigen::ColMajor> AtA(
       outer_product_->num_rows(),
       outer_product_->num_rows(),
@@ -186,14 +180,12 @@
 
   const Vector b = VectorRef(rhs_and_solution, outer_product_->num_rows());
   if (simplicial_ldlt_.get() == NULL || options_.dynamic_sparsity) {
-    typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double, Eigen::ColMajor>,
-                                  Eigen::Upper> SimplicialLDLT;
     simplicial_ldlt_.reset(new SimplicialLDLT);
     // This is a crappy way to be doing this. But right now Eigen does
     // not expose a way to do symbolic analysis with a given
     // permutation pattern, so we cannot use a block analysis of the
     // Jacobian.
-    simplicial_ldlt_->analyzePattern(AtA.selfadjointView<Eigen::Upper>());
+    simplicial_ldlt_->analyzePattern(AtA);
     if (simplicial_ldlt_->info() != Eigen::Success) {
       summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
       summary.message =
@@ -203,7 +195,7 @@
   }
   event_logger.AddEvent("Analysis");
 
-  simplicial_ldlt_->factorize(AtA.selfadjointView<Eigen::Upper>());
+  simplicial_ldlt_->factorize(AtA);
   if(simplicial_ldlt_->info() != Eigen::Success) {
     summary.termination_type = LINEAR_SOLVER_FAILURE;
     summary.message =

internal/ceres/sparse_normal_cholesky_solver.h

diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index 0e6a988..6572835 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -93,9 +93,9 @@
   cs_dis* cxsparse_factor_;
 
 #ifdef CERES_USE_EIGEN_SPARSE
-  scoped_ptr<Eigen::SimplicialLDLT<
-               Eigen::SparseMatrix<double, Eigen::ColMajor>,
-               Eigen::Upper> > simplicial_ldlt_;
+  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
+                                Eigen::Upper> SimplicialLDLT;
+  scoped_ptr<SimplicialLDLT> simplicial_ldlt_;
 #endif
 
   scoped_ptr<CompressedRowSparseMatrix> outer_product_;