Add SparseCholesky

SparseCholesky is an interface to sparse cholesky factorization
routines across sparse linear algebra libraries. Each sparse
linear algebra library is responsible for implementing its own
instance of this interface.

As a result the various places - SparseNormalCholeskySolver,
SparseSchurComplementSolver and VisibilityBasedPreconditioner
are significantly simplified.

Change-Id: I8b465705eae83bba9e1adfffcc741a05c70faf2e

internal/ceres/CMakeLists.txt

@@ -63,6 +63,7 @@ set(CERES_INTERNAL_SRC
     dynamic_compressed_row_sparse_matrix.cc
     dynamic_sparse_normal_cholesky_solver.cc
     evaluator.cc
+    eigensparse.cc
     file.cc
     gradient_checker.cc
     gradient_checking_cost_function.cc
@@ -105,6 +106,7 @@ set(CERES_INTERNAL_SRC
     solver.cc
     solver_utils.cc
     sparse_matrix.cc
+    sparse_cholesky.cc
     sparse_normal_cholesky_solver.cc
     split.cc
     stringprintf.cc
@@ -339,6 +341,7 @@ if (BUILD_TESTING AND GFLAGS)
   ceres_test(single_linkage_clustering)
   ceres_test(small_blas)
   ceres_test(solver)
+  ceres_test(sparse_cholesky)
 
   # TODO(sameeragarwal): This test should ultimately be made
   # independent of SuiteSparse.

internal/ceres/canonical_views_clustering.cc

@@ -29,11 +29,6 @@
 // Author: David Gallup (dgallup@google.com)
 //         Sameer Agarwal (sameeragarwal@google.com)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/canonical_views_clustering.h"
 
 #include "ceres/collections_port.h"
@@ -243,5 +238,3 @@ void CanonicalViewsClustering::ComputeClusterMembership(
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/canonical_views_clustering.h

@@ -41,11 +41,6 @@
 #ifndef CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_
 #define CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include <vector>
 
 #include "ceres/collections_port.h"
@@ -132,5 +127,4 @@ struct CanonicalViewsClusteringOptions {
 }  // namespace internal
 }  // namespace ceres
 
-#endif  // CERES_NO_SUITESPARSE
 #endif  // CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_

internal/ceres/canonical_views_clustering_test.cc

@@ -29,11 +29,6 @@
 // Author: Sameer Agarwal (sameeragarwal@google.com)
 //         David Gallup (dgallup@google.com)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/canonical_views_clustering.h"
 
 #include "ceres/collections_port.h"
@@ -146,5 +141,3 @@ TEST_F(CanonicalViewsTest, SimilarityPenaltyTest) {
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/cxsparse.cc

@@ -35,7 +35,9 @@
 
 #include "ceres/cxsparse.h"
 
+#include <string>
 #include <vector>
+
 #include "ceres/compressed_col_sparse_matrix_utils.h"
 #include "ceres/compressed_row_sparse_matrix.h"
 #include "ceres/triplet_sparse_matrix.h"
@@ -46,8 +48,7 @@ namespace internal {
 
 using std::vector;
 
-CXSparse::CXSparse() : scratch_(NULL), scratch_size_(0) {
-}
+CXSparse::CXSparse() : scratch_(NULL), scratch_size_(0) {}
 
 CXSparse::~CXSparse() {
   if (scratch_size_ > 0) {
@@ -55,47 +56,42 @@ CXSparse::~CXSparse() {
   }
 }
 
+csn* CXSparse::Cholesky(cs_di* A, cs_dis* symbolic_factor) {
+  return cs_di_chol(A, symbolic_factor);
+}
 
-bool CXSparse::SolveCholesky(cs_di* A,
-                             cs_dis* symbolic_factorization,
-                             double* b) {
+void CXSparse::Solve(cs_dis* symbolic_factor, csn* numeric_factor, double* b) {
   // Make sure we have enough scratch space available.
-  if (scratch_size_ < A->n) {
+  const int num_cols = numeric_factor->L->n;
+  if (scratch_size_ < num_cols) {
     if (scratch_size_ > 0) {
       cs_di_free(scratch_);
     }
     scratch_ =
-        reinterpret_cast<CS_ENTRY*>(cs_di_malloc(A->n, sizeof(CS_ENTRY)));
-    scratch_size_ = A->n;
-  }
-
-  // Solve using Cholesky factorization
-  csn* numeric_factorization = cs_di_chol(A, symbolic_factorization);
-  if (numeric_factorization == NULL) {
-    LOG(WARNING) << "Cholesky factorization failed.";
-    return false;
+        reinterpret_cast<CS_ENTRY*>(cs_di_malloc(num_cols, sizeof(CS_ENTRY)));
+    scratch_size_ = num_cols;
   }
 
-  // When the Cholesky factorization succeeded, these methods are
+  // When the Cholesky factor succeeded, these methods are
   // guaranteed to succeeded as well. In the comments below, "x"
   // refers to the scratch space.
   //
   // Set x = P * b.
-  cs_di_ipvec(symbolic_factorization->pinv, b, scratch_, A->n);
+  cs_di_ipvec(symbolic_factor->pinv, b, scratch_, num_cols);
   // Set x = L \ x.
-  cs_di_lsolve(numeric_factorization->L, scratch_);
+  cs_di_lsolve(numeric_factor->L, scratch_);
   // Set x = L' \ x.
-  cs_di_ltsolve(numeric_factorization->L, scratch_);
+  cs_di_ltsolve(numeric_factor->L, scratch_);
   // Set b = P' * x.
-  cs_di_pvec(symbolic_factorization->pinv, scratch_, b, A->n);
+  cs_di_pvec(symbolic_factor->pinv, scratch_, b, num_cols);
+}
 
-  // Free Cholesky factorization.
-  cs_di_nfree(numeric_factorization);
-  return true;
+bool CXSparse::SolveCholesky(cs_di* lhs, double* rhs_and_solution) {
+  return cs_cholsol(1, lhs, rhs_and_solution);
 }
 
 cs_dis* CXSparse::AnalyzeCholesky(cs_di* A) {
-  // order = 1 for Cholesky factorization.
+  // order = 1 for Cholesky factor.
   return cs_schol(1, A);
 }
 
@@ -112,16 +108,12 @@ cs_dis* CXSparse::BlockAnalyzeCholesky(cs_di* A,
 
   vector<int> block_rows;
   vector<int> block_cols;
-  CompressedColumnScalarMatrixToBlockMatrix(A->i,
-                                            A->p,
-                                            row_blocks,
-                                            col_blocks,
-                                            &block_rows,
-                                            &block_cols);
+  CompressedColumnScalarMatrixToBlockMatrix(
+      A->i, A->p, row_blocks, col_blocks, &block_rows, &block_cols);
   cs_di block_matrix;
   block_matrix.m = num_row_blocks;
   block_matrix.n = num_col_blocks;
-  block_matrix.nz  = -1;
+  block_matrix.nz = -1;
   block_matrix.nzmax = block_rows.size();
   block_matrix.p = &block_cols[0];
   block_matrix.i = &block_rows[0];
@@ -135,34 +127,30 @@ cs_dis* CXSparse::BlockAnalyzeCholesky(cs_di* A,
   vector<int> scalar_ordering;
   BlockOrderingToScalarOrdering(row_blocks, block_ordering, &scalar_ordering);
 
-  cs_dis* symbolic_factorization =
+  cs_dis* symbolic_factor =
       reinterpret_cast<cs_dis*>(cs_calloc(1, sizeof(cs_dis)));
-  symbolic_factorization->pinv = cs_pinv(&scalar_ordering[0], A->n);
-  cs* permuted_A = cs_symperm(A, symbolic_factorization->pinv, 0);
-
-  symbolic_factorization->parent = cs_etree(permuted_A, 0);
-  int* postordering = cs_post(symbolic_factorization->parent, A->n);
-  int* column_counts = cs_counts(permuted_A,
-                                 symbolic_factorization->parent,
-                                 postordering,
-                                 0);
+  symbolic_factor->pinv = cs_pinv(&scalar_ordering[0], A->n);
+  cs* permuted_A = cs_symperm(A, symbolic_factor->pinv, 0);
+
+  symbolic_factor->parent = cs_etree(permuted_A, 0);
+  int* postordering = cs_post(symbolic_factor->parent, A->n);
+  int* column_counts =
+      cs_counts(permuted_A, symbolic_factor->parent, postordering, 0);
   cs_free(postordering);
   cs_spfree(permuted_A);
 
-  symbolic_factorization->cp = (int*) cs_malloc(A->n+1, sizeof(int));
-  symbolic_factorization->lnz = cs_cumsum(symbolic_factorization->cp,
-                                          column_counts,
-                                          A->n);
-  symbolic_factorization->unz = symbolic_factorization->lnz;
+  symbolic_factor->cp = (int*)cs_malloc(A->n + 1, sizeof(int));
+  symbolic_factor->lnz = cs_cumsum(symbolic_factor->cp, column_counts, A->n);
+  symbolic_factor->unz = symbolic_factor->lnz;
 
   cs_free(column_counts);
 
-  if (symbolic_factorization->lnz < 0) {
-    cs_sfree(symbolic_factorization);
-    symbolic_factorization = NULL;
+  if (symbolic_factor->lnz < 0) {
+    cs_sfree(symbolic_factor);
+    symbolic_factor = NULL;
   }
 
-  return symbolic_factorization;
+  return symbolic_factor;
 }
 
 cs_di CXSparse::CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A) {
@@ -196,20 +184,97 @@ void CXSparse::ApproximateMinimumDegreeOrdering(cs_di* A, int* ordering) {
   cs_free(cs_ordering);
 }
 
-cs_di* CXSparse::TransposeMatrix(cs_di* A) {
-  return cs_di_transpose(A, 1);
-}
+cs_di* CXSparse::TransposeMatrix(cs_di* A) { return cs_di_transpose(A, 1); }
 
 cs_di* CXSparse::MatrixMatrixMultiply(cs_di* A, cs_di* B) {
   return cs_di_multiply(A, B);
 }
 
-void CXSparse::Free(cs_di* sparse_matrix) {
-  cs_di_spfree(sparse_matrix);
+void CXSparse::Free(cs_di* sparse_matrix) { cs_di_spfree(sparse_matrix); }
+
+void CXSparse::Free(cs_dis* symbolic_factor) { cs_di_sfree(symbolic_factor); }
+
+void CXSparse::Free(csn* numeric_factor) { cs_di_nfree(numeric_factor); }
+
+CXSparseCholesky* CXSparseCholesky::Create(const OrderingType ordering_type) {
+  return new CXSparseCholesky(ordering_type);
+}
+
+CompressedRowSparseMatrix::StorageType CXSparseCholesky::StorageType() const {
+  return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
 }
 
-void CXSparse::Free(cs_dis* symbolic_factorization) {
-  cs_di_sfree(symbolic_factorization);
+CXSparseCholesky::CXSparseCholesky(const OrderingType ordering_type)
+    : ordering_type_(ordering_type),
+      symbolic_factor_(NULL),
+      numeric_factor_(NULL) {}
+
+CXSparseCholesky::~CXSparseCholesky() {
+  FreeSymbolicFactorization();
+  FreeNumericFactorization();
+}
+
+LinearSolverTerminationType CXSparseCholesky::Factorize(
+    CompressedRowSparseMatrix* lhs, std::string* message) {
+  CHECK_EQ(lhs->storage_type(), StorageType());
+  if (lhs == NULL) {
+    *message = "Failure: Input lhs is NULL.";
+    return LINEAR_SOLVER_FATAL_ERROR;
+  }
+
+  cs_di cs_lhs = cs_.CreateSparseMatrixTransposeView(lhs);
+
+  if (symbolic_factor_ == NULL) {
+    if (ordering_type_ == NATURAL) {
+      symbolic_factor_ = cs_.AnalyzeCholeskyWithNaturalOrdering(&cs_lhs);
+    } else {
+      if (!lhs->col_blocks().empty() && !(lhs->row_blocks().empty())) {
+        symbolic_factor_ = cs_.BlockAnalyzeCholesky(
+            &cs_lhs, lhs->col_blocks(), lhs->row_blocks());
+      } else {
+        symbolic_factor_ = cs_.AnalyzeCholesky(&cs_lhs);
+      }
+    }
+
+    if (symbolic_factor_ == NULL) {
+      *message = "CXSparse Failure : Symbolic factorization failed.";
+      return LINEAR_SOLVER_FATAL_ERROR;
+    }
+  }
+
+  FreeNumericFactorization();
+  numeric_factor_ = cs_.Cholesky(&cs_lhs, symbolic_factor_);
+  if (numeric_factor_ == NULL) {
+    *message = "CXSparse Failure : Numeric factorization failed.";
+    return LINEAR_SOLVER_FAILURE;
+  }
+
+  return LINEAR_SOLVER_SUCCESS;
+}
+
+LinearSolverTerminationType CXSparseCholesky::Solve(const double* rhs,
+                                                    double* solution,
+                                                    std::string* message) {
+  CHECK(numeric_factor_ != NULL)
+      << "Solve called without a call to Factorize first.";
+  const int num_cols = numeric_factor_->L->n;
+  memcpy(solution, rhs, num_cols * sizeof(*solution));
+  cs_.Solve(symbolic_factor_, numeric_factor_, solution);
+  return LINEAR_SOLVER_SUCCESS;
+}
+
+void CXSparseCholesky::FreeSymbolicFactorization() {
+  if (symbolic_factor_ != NULL) {
+    cs_.Free(symbolic_factor_);
+    symbolic_factor_ = NULL;
+  }
+}
+
+void CXSparseCholesky::FreeNumericFactorization() {
+  if (numeric_factor_ != NULL) {
+    cs_.Free(numeric_factor_);
+    numeric_factor_ = NULL;
+  }
 }
 
 }  // namespace internal

internal/ceres/cxsparse.h

@@ -36,7 +36,11 @@
 
 #ifndef CERES_NO_CXSPARSE
 
+#include <string>
 #include <vector>
+
+#include "ceres/linear_solver.h"
+#include "ceres/sparse_cholesky.h"
 #include "cs.h"
 
 namespace ceres {
@@ -54,14 +58,18 @@ class CXSparse {
   CXSparse();
   ~CXSparse();
 
-  // Solves a symmetric linear system A * x = b using Cholesky factorization.
-  //  A                      - The system matrix.
-  //  symbolic_factorization - The symbolic factorization of A. This is obtained
-  //                           from AnalyzeCholesky.
-  //  b                      - The right hand size of the linear equation. This
-  //                           array will also recieve the solution.
-  // Returns false if Cholesky factorization of A fails.
-  bool SolveCholesky(cs_di* A, cs_dis* symbolic_factorization, double* b);
+  // Solve the system lhs * solution = rhs in place by using an
+  // approximate minimum degree fill reducing ordering.
+  bool SolveCholesky(cs_di* lhs, double* rhs_and_solution);
+
+  // Solves a linear system given its symbolic and numeric factorization.
+  void Solve(cs_dis* symbolic_factor, csn* numeric_factor, double* rhs_and_solution);
+
+  // Compute the numeric Cholesky factorization of A, given its
+  // symbolic factorization.
+  //
+  // Caller owns the result.
+  csn* Cholesky(cs_di* A, cs_dis* symbolic_factor);
 
   // Creates a sparse matrix from a compressed-column form. No memory is
   // allocated or copied; the structure A is filled out with info from the
@@ -117,6 +125,7 @@ class CXSparse {
 
   void Free(cs_di* sparse_matrix);
   void Free(cs_dis* symbolic_factorization);
+  void Free(csn* numeric_factorization);
 
  private:
   // Cached scratch space
@@ -124,6 +133,30 @@ class CXSparse {
   int scratch_size_;
 };
 
+class CXSparseCholesky : public SparseCholesky {
+ public:
+  // Factory
+  static CXSparseCholesky* Create(const OrderingType ordering_type);
+
+  // SparseCholesky interface.
+  virtual ~CXSparseCholesky();
+  virtual CompressedRowSparseMatrix::StorageType StorageType() const;
+  virtual LinearSolverTerminationType Factorize(
+      CompressedRowSparseMatrix* lhs, std::string* message);
+  virtual LinearSolverTerminationType Solve(const double* rhs,
+                                            double* solution,
+                                            std::string* message);
+ private:
+  CXSparseCholesky(const OrderingType ordering_type);
+  void FreeSymbolicFactorization();
+  void FreeNumericFactorization();
+
+  const OrderingType ordering_type_;
+  CXSparse cs_;
+  cs_dis* symbolic_factor_;
+  csn* numeric_factor_;
+};
+
 }  // namespace internal
 }  // namespace ceres
 

internal/ceres/dynamic_sparse_normal_cholesky_solver.cc

@@ -211,20 +211,10 @@ LinearSolver::Summary DynamicSparseNormalCholeskySolver::SolveImplUsingCXSparse(
   cxsparse.Free(a);
   event_logger.AddEvent("NormalEquations");
 
-  cs_dis* factor = cxsparse.AnalyzeCholesky(lhs);
-  event_logger.AddEvent("Analysis");
-
-  if (factor == NULL) {
-    summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-    summary.message = "CXSparse::AnalyzeCholesky failed.";
-  } else if (!cxsparse.SolveCholesky(lhs, factor, rhs_and_solution)) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "CXSparse::SolveCholesky failed.";
-  }
+  cxsparse.SolveCholesky(lhs, rhs_and_solution);
   event_logger.AddEvent("Solve");
 
   cxsparse.Free(lhs);
-  cxsparse.Free(factor);
   event_logger.AddEvent("TearDown");
   return summary;
 #endif

internal/ceres/eigensparse.cc

@@ -0,0 +1,143 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/eigensparse.h"
+
+#ifdef CERES_USE_EIGEN_SPARSE
+
+#include <sstream>
+#include "Eigen/SparseCholesky"
+#include "Eigen/SparseCore"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/linear_solver.h"
+
+namespace ceres {
+namespace internal {
+
+template <typename Solver>
+class EigenSparseCholeskyTemplate : public EigenSparseCholesky {
+ public:
+  EigenSparseCholeskyTemplate() : analyzed_(false) {}
+  virtual ~EigenSparseCholeskyTemplate() {}
+  virtual CompressedRowSparseMatrix::StorageType StorageType() const {
+    return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
+  }
+
+  virtual LinearSolverTerminationType Factorize(
+      const Eigen::SparseMatrix<double>& lhs, std::string* message) {
+    if (!analyzed_) {
+      solver_.analyzePattern(lhs);
+
+      if (VLOG_IS_ON(2)) {
+        std::stringstream ss;
+        solver_.dumpMemory(ss);
+        VLOG(2) << "Symbolic Analysis\n" << ss.str();
+      }
+
+      if (solver_.info() != Eigen::Success) {
+        *message = "Eigen failure. Unable to find symbolic factorization.";
+        return LINEAR_SOLVER_FATAL_ERROR;
+      }
+
+      analyzed_ = true;
+    }
+
+    solver_.factorize(lhs);
+    if (solver_.info() != Eigen::Success) {
+      *message = "Eigen failure. Unable to find numeric factorization.";
+      return LINEAR_SOLVER_FAILURE;
+    }
+    return LINEAR_SOLVER_SUCCESS;
+  }
+
+  virtual LinearSolverTerminationType Solve(const double* rhs,
+                                            double* solution,
+                                            std::string* message) {
+    CHECK(analyzed_) << "Solve called without a call to Factorize first.";
+
+    VectorRef(solution, solver_.cols()) =
+        solver_.solve(ConstVectorRef(rhs, solver_.cols()));
+    if (solver_.info() != Eigen::Success) {
+      *message = "Eigen failure. Unable to do triangular solve.";
+      return LINEAR_SOLVER_FAILURE;
+    }
+    return LINEAR_SOLVER_SUCCESS;
+  }
+
+  virtual LinearSolverTerminationType Factorize(CompressedRowSparseMatrix* lhs,
+                                                std::string* message) {
+    CHECK_EQ(lhs->storage_type(), StorageType());
+    Eigen::MappedSparseMatrix<double, Eigen::ColMajor> eigen_lhs(
+        lhs->num_rows(),
+        lhs->num_rows(),
+        lhs->num_nonzeros(),
+        lhs->mutable_rows(),
+        lhs->mutable_cols(),
+        lhs->mutable_values());
+    return Factorize(eigen_lhs, message);
+  }
+
+ private:
+  bool analyzed_;
+  Solver solver_;
+};
+
+EigenSparseCholesky* EigenSparseCholesky::Create(
+    const OrderingType ordering_type) {
+  // The preprocessor gymnastics here are dealing with the fact that
+  // before version 3.2.2, Eigen did not support a third template
+  // parameter to specify the ordering and it always defaults to AMD.
+#if EIGEN_VERSION_AT_LEAST(3, 2, 2)
+  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
+                                Eigen::Upper,
+                                Eigen::AMDOrdering<int> >
+      WithAMDOrdering;
+  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
+                                Eigen::Upper,
+                                Eigen::NaturalOrdering<int> >
+      WithNaturalOrdering;
+  if (ordering_type == AMD) {
+    return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
+  } else {
+    return new EigenSparseCholeskyTemplate<WithNaturalOrdering>();
+  }
+#else
+  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper>
+      WithAMDOrdering;
+  return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
+#endif
+}
+
+EigenSparseCholesky::~EigenSparseCholesky() {}
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_USE_EIGEN_SPARSE

internal/ceres/eigensparse.h

@@ -0,0 +1,69 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// A simple C++ interface to the Eigen's Sparse Cholesky routines.
+
+#ifndef CERES_INTERNAL_EIGENSPARSE_H_
+#define CERES_INTERNAL_EIGENSPARSE_H_
+
+// This include must come before any #ifndef check on Ceres compile options.
+#include "ceres/internal/port.h"
+
+#ifdef CERES_USE_EIGEN_SPARSE
+
+#include <string>
+#include "Eigen/SparseCore"
+#include "ceres/linear_solver.h"
+#include "ceres/sparse_cholesky.h"
+
+namespace ceres {
+namespace internal {
+
+class EigenSparseCholesky : public SparseCholesky {
+ public:
+  // Factory
+  static EigenSparseCholesky* Create(const OrderingType ordering_type);
+
+  // SparseCholesky interface.
+  virtual ~EigenSparseCholesky();
+  virtual CompressedRowSparseMatrix::StorageType StorageType() const = 0;
+  virtual LinearSolverTerminationType Factorize(
+      const Eigen::SparseMatrix<double>& lhs, std::string* message) = 0;
+  virtual LinearSolverTerminationType Solve(const double* rhs,
+                                            double* solution,
+                                            std::string* message) = 0;
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_USE_EIGEN_SPARSE
+
+#endif  // CERES_INTERNAL_EIGENSPARSE_H_

internal/ceres/linear_solver.h

@@ -69,6 +69,10 @@ enum LinearSolverTerminationType {
   LINEAR_SOLVER_FATAL_ERROR
 };
 
+enum OrderingType {
+  NATURAL,
+  AMD
+};
 
 class LinearOperator;
 

internal/ceres/schur_complement_solver.cc

@@ -28,33 +28,30 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
-#include "ceres/internal/port.h"
+#include "ceres/schur_complement_solver.h"
 
 #include <algorithm>
 #include <ctime>
 #include <set>
-#include <sstream>
 #include <vector>
 
+#include "Eigen/Dense"
+#include "Eigen/SparseCore"
 #include "ceres/block_random_access_dense_matrix.h"
 #include "ceres/block_random_access_matrix.h"
 #include "ceres/block_random_access_sparse_matrix.h"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
 #include "ceres/conjugate_gradients_solver.h"
-#include "ceres/cxsparse.h"
 #include "ceres/detect_structure.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/lapack.h"
 #include "ceres/linear_solver.h"
-#include "ceres/schur_complement_solver.h"
-#include "ceres/suitesparse.h"
+#include "ceres/sparse_cholesky.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "ceres/types.h"
 #include "ceres/wall_time.h"
-#include "Eigen/Dense"
-#include "Eigen/SparseCore"
 
 namespace ceres {
 namespace internal {
@@ -139,6 +136,7 @@ LinearSolver::Summary SchurComplementSolver::SolveImpl(
     eliminator_->Init(
         options_.elimination_groups[0], kFullRankETE, A->block_structure());
   };
+
   std::fill(x, x + A->num_cols(), 0.0);
   event_logger.AddEvent("Setup");
 
@@ -227,21 +225,13 @@ DenseSchurComplementSolver::SolveReducedLinearSystem(
 
 SparseSchurComplementSolver::SparseSchurComplementSolver(
     const LinearSolver::Options& options)
-    : SchurComplementSolver(options),
-      factor_(NULL),
-      cxsparse_factor_(NULL) {
+    : SchurComplementSolver(options) {
+  sparse_cholesky_.reset(
+      SparseCholesky::Create(options.sparse_linear_algebra_library_type,
+                             options.use_postordering ? AMD : NATURAL));
 }
 
 SparseSchurComplementSolver::~SparseSchurComplementSolver() {
-  if (factor_ != NULL) {
-    ss_.Free(factor_);
-    factor_ = NULL;
-  }
-
-  if (cxsparse_factor_ != NULL) {
-    cxsparse_.Free(cxsparse_factor_);
-    cxsparse_factor_ = NULL;
-  }
 }
 
 // Determine the non-zero blocks in the Schur Complement matrix, and
@@ -315,297 +305,47 @@ void SparseSchurComplementSolver::InitStorage(
   set_rhs(new double[lhs()->num_rows()]);
 }
 
-LinearSolver::Summary
-SparseSchurComplementSolver::SolveReducedLinearSystem(
-          const LinearSolver::PerSolveOptions& per_solve_options,
-          double* solution) {
+LinearSolver::Summary SparseSchurComplementSolver::SolveReducedLinearSystem(
+    const LinearSolver::PerSolveOptions& per_solve_options, double* solution) {
   if (options().type == ITERATIVE_SCHUR) {
-    CHECK(options().use_explicit_schur_complement);
     return SolveReducedLinearSystemUsingConjugateGradients(per_solve_options,
                                                            solution);
   }
 
-  switch (options().sparse_linear_algebra_library_type) {
-    case SUITE_SPARSE:
-      return SolveReducedLinearSystemUsingSuiteSparse(per_solve_options,
-                                                      solution);
-    case CX_SPARSE:
-      return SolveReducedLinearSystemUsingCXSparse(per_solve_options,
-                                                   solution);
-    case EIGEN_SPARSE:
-      return SolveReducedLinearSystemUsingEigen(per_solve_options,
-                                                solution);
-    default:
-      LOG(FATAL) << "Unknown sparse linear algebra library : "
-                 << options().sparse_linear_algebra_library_type;
-  }
-
-  return LinearSolver::Summary();
-}
-
-// Solve the system Sx = r, assuming that the matrix S is stored in a
-// BlockRandomAccessSparseMatrix.  The linear system is solved using
-// CHOLMOD's sparse cholesky factorization routines.
-LinearSolver::Summary
-SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
-    const LinearSolver::PerSolveOptions& per_solve_options,
-    double* solution) {
-#ifdef CERES_NO_SUITESPARSE
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-  summary.message = "Ceres was not built with SuiteSparse support. "
-      "Therefore, SPARSE_SCHUR cannot be used with SUITE_SPARSE";
-  return summary;
-
-#else
-
   LinearSolver::Summary summary;
   summary.num_iterations = 0;
   summary.termination_type = LINEAR_SOLVER_SUCCESS;
   summary.message = "Success.";
 
-  TripletSparseMatrix* tsm =
-      const_cast<TripletSparseMatrix*>(
-          down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
-  const int num_rows = tsm->num_rows();
-
-  // The case where there are no f blocks, and the system is block
-  // diagonal.
-  if (num_rows == 0) {
+  const TripletSparseMatrix* tsm =
+      down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix();
+  if (tsm->num_rows() == 0) {
     return summary;
   }
 
-  summary.num_iterations = 1;
-  cholmod_sparse* cholmod_lhs = NULL;
-  if (options().use_postordering) {
-    // If we are going to do a full symbolic analysis of the schur
-    // complement matrix from scratch and not rely on the
-    // pre-ordering, then the fastest path in cholmod_factorize is the
-    // one corresponding to upper triangular matrices.
-
-    // Create a upper triangular symmetric matrix.
-    cholmod_lhs = ss_.CreateSparseMatrix(tsm);
-    cholmod_lhs->stype = 1;
-
-    if (factor_ == NULL) {
-      factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs,
-                                         blocks_,
-                                         blocks_,
-                                         &summary.message);
-    }
+  scoped_ptr<CompressedRowSparseMatrix> lhs;
+  const CompressedRowSparseMatrix::StorageType storage_type =
+      sparse_cholesky_->StorageType();
+  if (storage_type == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+    lhs.reset(CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
+    lhs->set_storage_type(CompressedRowSparseMatrix::UPPER_TRIANGULAR);
   } else {
-    // If we are going to use the natural ordering (i.e. rely on the
-    // pre-ordering computed by solver_impl.cc), then the fastest
-    // path in cholmod_factorize is the one corresponding to lower
-    // triangular matrices.
-
-    // Create a upper triangular symmetric matrix.
-    cholmod_lhs = ss_.CreateSparseMatrixTranspose(tsm);
-    cholmod_lhs->stype = -1;
-
-    if (factor_ == NULL) {
-      factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(cholmod_lhs,
-                                                       &summary.message);
-    }
-  }
-
-  if (factor_ == NULL) {
-    ss_.Free(cholmod_lhs);
-    summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-    // No need to set message as it has already been set by the
-    // symbolic analysis routines above.
-    return summary;
-  }
-
-  summary.termination_type =
-    ss_.Cholesky(cholmod_lhs, factor_, &summary.message);
-
-  ss_.Free(cholmod_lhs);
-
-  if (summary.termination_type != LINEAR_SOLVER_SUCCESS) {
-    // No need to set message as it has already been set by the
-    // numeric factorization routine above.
-    return summary;
-  }
-
-  cholmod_dense*  cholmod_rhs =
-      ss_.CreateDenseVector(const_cast<double*>(rhs()), num_rows, num_rows);
-  cholmod_dense* cholmod_solution = ss_.Solve(factor_,
-                                              cholmod_rhs,
-                                              &summary.message);
-  ss_.Free(cholmod_rhs);
-
-  if (cholmod_solution == NULL) {
-    summary.message =
-        "SuiteSparse failure. Unable to perform triangular solve.";
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    return summary;
-  }
-
-  VectorRef(solution, num_rows)
-      = VectorRef(static_cast<double*>(cholmod_solution->x), num_rows);
-  ss_.Free(cholmod_solution);
-  return summary;
-#endif  // CERES_NO_SUITESPARSE
-}
-
-// Solve the system Sx = r, assuming that the matrix S is stored in a
-// BlockRandomAccessSparseMatrix.  The linear system is solved using
-// CXSparse's sparse cholesky factorization routines.
-LinearSolver::Summary
-SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse(
-    const LinearSolver::PerSolveOptions& per_solve_options,
-    double* solution) {
-#ifdef CERES_NO_CXSPARSE
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-  summary.message = "Ceres was not built with CXSparse support. "
-      "Therefore, SPARSE_SCHUR cannot be used with CX_SPARSE";
-  return summary;
-
-#else
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  summary.message = "Success.";
-
-  // Extract the TripletSparseMatrix that is used for actually storing S.
-  TripletSparseMatrix* tsm =
-      const_cast<TripletSparseMatrix*>(
-          down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
-  const int num_rows = tsm->num_rows();
-
-  // The case where there are no f blocks, and the system is block
-  // diagonal.
-  if (num_rows == 0) {
-    return summary;
-  }
-
-  cs_di* lhs = CHECK_NOTNULL(cxsparse_.CreateSparseMatrix(tsm));
-  VectorRef(solution, num_rows) = ConstVectorRef(rhs(), num_rows);
-
-  // Compute symbolic factorization if not available.
-  if (cxsparse_factor_ == NULL) {
-    cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(lhs, blocks_, blocks_);
-  }
-
-  if (cxsparse_factor_ == NULL) {
-    summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-    summary.message =
-        "CXSparse failure. Unable to find symbolic factorization.";
-  } else if (!cxsparse_.SolveCholesky(lhs, cxsparse_factor_, solution)) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "CXSparse::SolveCholesky failed.";
-  }
-
-  cxsparse_.Free(lhs);
-  return summary;
-#endif  // CERES_NO_CXPARSE
-}
-
-// Solve the system Sx = r, assuming that the matrix S is stored in a
-// BlockRandomAccessSparseMatrix.  The linear system is solved using
-// Eigen's sparse cholesky factorization routines.
-LinearSolver::Summary
-SparseSchurComplementSolver::SolveReducedLinearSystemUsingEigen(
-    const LinearSolver::PerSolveOptions& per_solve_options,
-    double* solution) {
-#ifndef CERES_USE_EIGEN_SPARSE
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-  summary.message =
-      "SPARSE_SCHUR cannot be used with EIGEN_SPARSE. "
-      "Ceres was not built with support for "
-      "Eigen's SimplicialLDLT decomposition. "
-      "This requires enabling building with -DEIGENSPARSE=ON.";
-  return summary;
-
-#else
-  EventLogger event_logger("SchurComplementSolver::EigenSolve");
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  summary.message = "Success.";
-
-  // Extract the TripletSparseMatrix that is used for actually storing S.
-  TripletSparseMatrix* tsm =
-      const_cast<TripletSparseMatrix*>(
-          down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
-  const int num_rows = tsm->num_rows();
-
-  // The case where there are no f blocks, and the system is block
-  // diagonal.
-  if (num_rows == 0) {
-    return summary;
-  }
-
-  // This is an upper triangular matrix.
-  scoped_ptr<CompressedRowSparseMatrix> crsm(
-      CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
-  // Map this to a column major, lower triangular matrix.
-  Eigen::MappedSparseMatrix<double, Eigen::ColMajor> eigen_lhs(
-      crsm->num_rows(),
-      crsm->num_rows(),
-      crsm->num_nonzeros(),
-      crsm->mutable_rows(),
-      crsm->mutable_cols(),
-      crsm->mutable_values());
-  event_logger.AddEvent("ToCompressedRowSparseMatrix");
-
-  // Compute symbolic factorization if one does not exist.
-  if (simplicial_ldlt_.get() == NULL) {
-    simplicial_ldlt_.reset(new SimplicialLDLT);
-    // This ordering is quite bad. The scalar ordering produced by the
-    // AMD algorithm is quite bad and can be an order of magnitude
-    // worse than the one computed using the block version of the
-    // algorithm.
-    simplicial_ldlt_->analyzePattern(eigen_lhs);
-    if (VLOG_IS_ON(2)) {
-      std::stringstream ss;
-      simplicial_ldlt_->dumpMemory(ss);
-      VLOG(2) << "Symbolic Analysis\n"
-              << ss.str();
-    }
-    event_logger.AddEvent("Analysis");
-    if (simplicial_ldlt_->info() != Eigen::Success) {
-      summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-      summary.message =
-          "Eigen failure. Unable to find symbolic factorization.";
-      return summary;
-    }
-  }
-
-  simplicial_ldlt_->factorize(eigen_lhs);
-  event_logger.AddEvent("Factorize");
-  if (simplicial_ldlt_->info() != Eigen::Success) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "Eigen failure. Unable to find numeric factoriztion.";
-    return summary;
-  }
-
-  VectorRef(solution, num_rows) =
-      simplicial_ldlt_->solve(ConstVectorRef(rhs(), num_rows));
-  event_logger.AddEvent("Solve");
-  if (simplicial_ldlt_->info() != Eigen::Success) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "Eigen failure. Unable to do triangular solve.";
+    lhs.reset(
+        CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm));
+    lhs->set_storage_type(CompressedRowSparseMatrix::LOWER_TRIANGULAR);
   }
 
+  summary.num_iterations = 1;
+  summary.termination_type = sparse_cholesky_->FactorAndSolve(
+      lhs.get(), rhs(), solution, &summary.message);
   return summary;
-#endif  // CERES_USE_EIGEN_SPARSE
 }
 
 LinearSolver::Summary
 SparseSchurComplementSolver::SolveReducedLinearSystemUsingConjugateGradients(
     const LinearSolver::PerSolveOptions& per_solve_options,
     double* solution) {
+  CHECK(options().use_explicit_schur_complement);
   const int num_rows = lhs()->num_rows();
   // The case where there are no f blocks, and the system is block
   // diagonal.

internal/ceres/schur_complement_solver.h

@@ -35,18 +35,15 @@
 #include <utility>
 #include <vector>
 
-#include "ceres/internal/port.h"
-
+#include "ceres/block_random_access_diagonal_matrix.h"
 #include "ceres/block_random_access_matrix.h"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
-#include "ceres/cxsparse.h"
+#include "ceres/internal/port.h"
+#include "ceres/internal/scoped_ptr.h"
 #include "ceres/linear_solver.h"
 #include "ceres/schur_eliminator.h"
-#include "ceres/suitesparse.h"
-#include "ceres/internal/scoped_ptr.h"
 #include "ceres/types.h"
-#include "ceres/block_random_access_diagonal_matrix.h"
 
 #ifdef CERES_USE_EIGEN_SPARSE
 #include "Eigen/SparseCholesky"
@@ -57,6 +54,7 @@ namespace ceres {
 namespace internal {
 
 class BlockSparseMatrix;
+class SparseCholesky;
 
 // Base class for Schur complement based linear least squares
 // solvers. It assumes that the input linear system Ax = b can be
@@ -96,18 +94,19 @@ class BlockSparseMatrix;
 // be used for problems with upto a few hundred cameras.
 //
 // SparseSchurComplementSolver: For problems where the Schur
-// complement matrix is large and sparse. It requires that
-// CHOLMOD/SuiteSparse be installed, as it uses CHOLMOD to find a
-// sparse Cholesky factorization of the Schur complement. This solver
-// can be used for solving structure from motion problems with tens of
-// thousands of cameras, though depending on the exact sparsity
-// structure, it maybe better to use an iterative solver.
+// complement matrix is large and sparse. It requires that Ceres be
+// build with at least one sparse linear algebra library, as it
+// computes a sparse Cholesky factorization of the Schur complement.
+//
+// This solver can be used for solving structure from motion problems
+// with tens of thousands of cameras, though depending on the exact
+// sparsity structure, it maybe better to use an iterative solver.
 //
 // The two solvers can be instantiated by calling
 // LinearSolver::CreateLinearSolver with LinearSolver::Options::type
 // set to DENSE_SCHUR and SPARSE_SCHUR
-// respectively. LinearSolver::Options::elimination_groups[0] should be
-// at least 1.
+// respectively. LinearSolver::Options::elimination_groups[0] should
+// be at least 1.
 class SchurComplementSolver : public BlockSparseMatrixSolver {
  public:
   explicit SchurComplementSolver(const LinearSolver::Options& options)
@@ -174,48 +173,13 @@ class SparseSchurComplementSolver : public SchurComplementSolver {
   virtual LinearSolver::Summary SolveReducedLinearSystem(
       const LinearSolver::PerSolveOptions& per_solve_options,
       double* solution);
-  LinearSolver::Summary SolveReducedLinearSystemUsingSuiteSparse(
-      const LinearSolver::PerSolveOptions& per_solve_options,
-      double* solution);
-  LinearSolver::Summary SolveReducedLinearSystemUsingCXSparse(
-      const LinearSolver::PerSolveOptions& per_solve_options,
-      double* solution);
-  LinearSolver::Summary SolveReducedLinearSystemUsingEigen(
-      const LinearSolver::PerSolveOptions& per_solve_options,
-      double* solution);
   LinearSolver::Summary SolveReducedLinearSystemUsingConjugateGradients(
       const LinearSolver::PerSolveOptions& per_solve_options,
       double* solution);
 
   // Size of the blocks in the Schur complement.
   std::vector<int> blocks_;
-
-  SuiteSparse ss_;
-  // Symbolic factorization of the reduced linear system. Precomputed
-  // once and reused in subsequent calls.
-  cholmod_factor* factor_;
-
-  CXSparse cxsparse_;
-  // Cached factorization
-  cs_dis* cxsparse_factor_;
-
-#ifdef CERES_USE_EIGEN_SPARSE
-
-  // The preprocessor gymnastics here are dealing with the fact that
-  // before version 3.2.2, Eigen did not support a third template
-  // parameter to specify the ordering.
-#if EIGEN_VERSION_AT_LEAST(3,2,2)
-  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Lower,
-                                Eigen::NaturalOrdering<int> >
-  SimplicialLDLT;
-#else
-  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Lower>
-  SimplicialLDLT;
-#endif
-
-  scoped_ptr<SimplicialLDLT> simplicial_ldlt_;
-#endif
-
+  scoped_ptr<SparseCholesky> sparse_cholesky_;
   scoped_ptr<BlockRandomAccessDiagonalMatrix> preconditioner_;
   CERES_DISALLOW_COPY_AND_ASSIGN(SparseSchurComplementSolver);
 };

internal/ceres/schur_complement_solver_test.cc

@@ -28,7 +28,10 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
+#include "ceres/schur_complement_solver.h"
+
 #include <cstddef>
+
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
 #include "ceres/casts.h"
@@ -36,7 +39,6 @@
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/linear_least_squares_problems.h"
 #include "ceres/linear_solver.h"
-#include "ceres/schur_complement_solver.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "ceres/types.h"
 #include "glog/logging.h"
@@ -120,12 +122,13 @@ class SchurComplementSolverTest : public ::testing::Test {
     EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
 
     if (regularization) {
+
       ASSERT_NEAR((sol_d - x).norm() / num_cols, 0, 1e-10)
-          << "Expected solution: " << sol_d.transpose()
+          << "Regularized Expected solution: " << sol_d.transpose()
           << " Actual solution: " << x.transpose();
     } else {
       ASSERT_NEAR((sol - x).norm() / num_cols, 0, 1e-10)
-          << "Expected solution: " << sol.transpose()
+          << "Unregularized Expected solution: " << sol.transpose()
           << " Actual solution: " << x.transpose();
     }
   }
@@ -142,27 +145,29 @@ class SchurComplementSolverTest : public ::testing::Test {
   Vector sol_d;
 };
 
-TEST_F(SchurComplementSolverTest, EigenBasedDenseSchurWithSmallProblem) {
+// TODO(sameeragarwal): Refactor these using value parameterized tests.
+// TODO(sameeragarwal): More extensive tests using random matrices.
+TEST_F(SchurComplementSolverTest, DenseSchurWithEigenSmallProblem) {
   ComputeAndCompareSolutions(2, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
   ComputeAndCompareSolutions(2, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 }
 
-TEST_F(SchurComplementSolverTest, EigenBasedDenseSchurWithLargeProblem) {
+TEST_F(SchurComplementSolverTest, DenseSchurWithEigenLargeProblem) {
   ComputeAndCompareSolutions(3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
   ComputeAndCompareSolutions(3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 }
 
-TEST_F(SchurComplementSolverTest, EigenBasedDenseSchurWithVaryingFBlockSize) {
+TEST_F(SchurComplementSolverTest, DenseSchurWithEigenVaryingFBlockSize) {
   ComputeAndCompareSolutions(4, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 }
 
 #ifndef CERES_NO_LAPACK
-TEST_F(SchurComplementSolverTest, LAPACKBasedDenseSchurWithSmallProblem) {
+TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKSmallProblem) {
   ComputeAndCompareSolutions(2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
   ComputeAndCompareSolutions(2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
 }
 
-TEST_F(SchurComplementSolverTest, LAPACKBasedDenseSchurWithLargeProblem) {
+TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKLargeProblem) {
   ComputeAndCompareSolutions(3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
   ComputeAndCompareSolutions(3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
 }

internal/ceres/single_linkage_clustering.cc

@@ -28,11 +28,6 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/single_linkage_clustering.h"
 
 #include "ceres/graph.h"
@@ -106,5 +101,3 @@ int ComputeSingleLinkageClustering(
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/single_linkage_clustering.h

@@ -31,11 +31,6 @@
 #ifndef CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_
 #define CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/collections_port.h"
 #include "ceres/graph.h"
 
@@ -70,5 +65,4 @@ int ComputeSingleLinkageClustering(
 }  // namespace internal
 }  // namespace ceres
 
-#endif  // CERES_NO_SUITESPARSE
 #endif  // CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_

internal/ceres/single_linkage_clustering_test.cc

@@ -28,11 +28,6 @@
 //
 // Author: Sameer Agarwal (sameeragarwal@google.com)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/single_linkage_clustering.h"
 
 #include "ceres/collections_port.h"
@@ -128,5 +123,3 @@ TEST(SingleLinkageClustering, ComponentWithWeakLinkAndStrongLink) {
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/sparse_cholesky.cc

@@ -0,0 +1,100 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/sparse_cholesky.h"
+
+#include "ceres/cxsparse.h"
+#include "ceres/eigensparse.h"
+#include "ceres/suitesparse.h"
+
+namespace ceres {
+namespace internal {
+
+SparseCholesky* SparseCholesky::Create(
+    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    OrderingType ordering_type) {
+  switch (sparse_linear_algebra_library_type) {
+    case SUITE_SPARSE:
+#ifndef CERES_NO_SUITESPARSE
+      return SuiteSparseCholesky::Create(ordering_type);
+#else
+      LOG(FATAL) << "Ceres was compiled without support for SuiteSparse.";
+      return NULL;
+#endif
+
+    case EIGEN_SPARSE:
+#ifdef CERES_USE_EIGEN_SPARSE
+      return EigenSparseCholesky::Create(ordering_type);
+#else
+      LOG(FATAL) << "Ceres was compiled without support for "
+                 << "Eigen's sparse Cholesky factorization routines.";
+      return NULL;
+#endif
+
+    case CX_SPARSE:
+#ifndef CERES_NO_CXSPARSE
+      return CXSparseCholesky::Create(ordering_type);
+#else
+      LOG(FATAL) << "Ceres was compiled without support for CXSparse.";
+      return NULL;
+#endif
+
+    default:
+      LOG(FATAL) << "Unknown sparse linear algebra library type : "
+                 << SparseLinearAlgebraLibraryTypeToString(
+                        sparse_linear_algebra_library_type);
+  }
+  return NULL;
+}
+
+SparseCholesky::~SparseCholesky() {}
+
+LinearSolverTerminationType SparseCholesky::FactorAndSolve(
+    CompressedRowSparseMatrix* lhs,
+    const double* rhs,
+    double* solution,
+    std::string* message) {
+  LinearSolverTerminationType termination_type = Factorize(lhs, message);
+  if (termination_type == LINEAR_SOLVER_SUCCESS) {
+    termination_type = Solve(rhs, solution, message);
+  }
+  return termination_type;
+}
+
+CompressedRowSparseMatrix::StorageType StorageTypeForSparseLinearAlgebraLibrary(
+    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type) {
+  if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
+    return CompressedRowSparseMatrix::UPPER_TRIANGULAR;
+  }
+  return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
+}
+
+}  // namespace internal
+}  // namespace ceres

internal/ceres/sparse_cholesky.h

@@ -0,0 +1,122 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_SPARSE_CHOLESKY_H_
+#define CERES_INTERNAL_SPARSE_CHOLESKY_H_
+
+// This include must come before any #ifndef check on Ceres compile options.
+#include "ceres/internal/port.h"
+
+#include "ceres/linear_solver.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+// An interface abstracts away the internal details of various sparse
+// linear algebra libraries and offers a simple API for solving
+// symmetric positive definite linear systems using a sparse Cholesky
+// factorization.
+//
+// Instances of SparseCholesky are expected to cache the symbolic
+// factorization of the linear system. They do this on the first call
+// to Factorize of FactorAndSolve. Subsequent calls to Factorize and
+// FactorAndSolve are expected to have the same sparsity structure.
+//
+// Example usage:
+//
+//  scoped_ptr<SparseCholesky>
+//  sparse_cholesky(SparseCholesky::Create(SUITE_SPARSE, AMD));
+//
+//  CompressedRowSparseMatrix lhs = ...;
+//  std::string message;
+//  CHECK_EQ(sparse_cholesky->Factorize(&lhs, &message), LINEAR_SOLVER_SUCCESS);
+//  Vector rhs = ...;
+//  Vector solution = ...;
+//  CHECK_EQ(sparse_cholesky->Solve(rhs.data(), solution.data(), &message),
+//           LINEAR_SOLVER_SUCCESS);
+
+class SparseCholesky {
+ public:
+  // Factory which returns an instance of SparseCholesky for the given
+  // sparse linear algebra library and fill reducing ordering
+  // strategy.
+  //
+  // Caller owns the result.
+  static SparseCholesky* Create(
+      SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+      OrderingType ordering_type);
+
+  virtual ~SparseCholesky();
+
+  // Due to the symmetry of the linear system, sparse linear algebra
+  // libraries only use one half of the input matrix. Whether it is
+  // the upper or the lower triangular part of the matrix depends on
+  // the library and the re-ordering strategy being used. This
+  // function tells the user the storage type expected of the input
+  // matrix for the sparse linear algebra library and reordering
+  // strategy used.
+  virtual CompressedRowSparseMatrix::StorageType StorageType() const = 0;
+
+  // Compute the numeric factorization of the given matrix.  If this
+  // is the first call to Factorize, first the symbolic factorization
+  // will be computed and cached and the numeric factorization will be
+  // computed based on that.
+  //
+  // Subsequent calls to Factorize will use that symbolic
+  // factorization assuming that the sparsity of the matrix has
+  // remained constant.
+  virtual LinearSolverTerminationType Factorize(
+      CompressedRowSparseMatrix* lhs, std::string* message) = 0;
+
+  // Compute the solution to the equation
+  //
+  // lhs * solution = rhs
+  //
+  // rhs and solution can point to the same memory location.
+  virtual LinearSolverTerminationType Solve(const double* rhs,
+                                            double* solution,
+                                            std::string* message) = 0;
+
+  // Convenience method which combines a call to Factorize and
+  // Solve. Solve is only called is Factorize returns
+  // LINEAR_SOLVER_SUCCESS.
+  virtual LinearSolverTerminationType FactorAndSolve(
+      CompressedRowSparseMatrix* lhs,
+      const double* rhs,
+      double* solution,
+      std::string* message);
+
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_SPARSE_CHOLESKY_H_

internal/ceres/sparse_cholesky_test.cc

@@ -0,0 +1,217 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/sparse_cholesky.h"
+
+#include <numeric>
+#include <vector>
+
+#include "Eigen/Dense"
+#include "Eigen/SparseCore"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/random.h"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+CompressedRowSparseMatrix* CreateRandomSymmetricPositiveDefiniteMatrix(
+    const int num_col_blocks,
+    const int min_col_block_size,
+    const int max_col_block_size,
+    const double block_density,
+    const CompressedRowSparseMatrix::StorageType storage_type) {
+  // Create a random matrix
+  CompressedRowSparseMatrix::RandomMatrixOptions options;
+  options.num_col_blocks = num_col_blocks;
+  options.min_col_block_size = min_col_block_size;
+  options.max_col_block_size = max_col_block_size;
+
+  options.num_row_blocks = 2 * num_col_blocks;
+  options.min_row_block_size = 1;
+  options.max_row_block_size = max_col_block_size;
+  options.block_density = block_density;
+  scoped_ptr<CompressedRowSparseMatrix> random_crsm(
+      CompressedRowSparseMatrix::CreateRandomMatrix(options));
+
+  // Add a diagonal block sparse matrix to make it full rank.
+  Vector diagonal = Vector::Ones(random_crsm->num_cols());
+  scoped_ptr<CompressedRowSparseMatrix> block_diagonal(
+      CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
+          diagonal.data(), random_crsm->col_blocks()));
+  random_crsm->AppendRows(*block_diagonal);
+
+  // Compute output = random_crsm' * random_crsm
+  std::vector<int> program;
+  CompressedRowSparseMatrix* output =
+      CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
+          *random_crsm, storage_type, &program);
+  CompressedRowSparseMatrix::ComputeOuterProduct(*random_crsm, program, output);
+  return output;
+}
+
+bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
+                             const Vector& rhs,
+                             Vector* solution) {
+  Matrix eigen_lhs;
+  lhs.ToDenseMatrix(&eigen_lhs);
+  if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+    Matrix full_lhs = eigen_lhs.selfadjointView<Eigen::Upper>();
+    Eigen::LLT<Matrix, Eigen::Upper> llt =
+        eigen_lhs.selfadjointView<Eigen::Upper>().llt();
+    if (llt.info() != Eigen::Success) {
+      return false;
+    }
+    *solution = llt.solve(rhs);
+    return (llt.info() == Eigen::Success);
+  }
+
+  Matrix full_lhs = eigen_lhs.selfadjointView<Eigen::Lower>();
+  Eigen::LLT<Matrix, Eigen::Lower> llt =
+      eigen_lhs.selfadjointView<Eigen::Lower>().llt();
+  if (llt.info() != Eigen::Success) {
+    return false;
+  }
+  *solution = llt.solve(rhs);
+  return (llt.info() == Eigen::Success);
+}
+
+void SparseCholeskySolverUnitTest(
+    const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const OrderingType ordering_type,
+    const bool use_block_structure,
+    const int num_blocks,
+    const int min_block_size,
+    const int max_block_size,
+    const double block_density) {
+  scoped_ptr<SparseCholesky> sparse_cholesky(SparseCholesky::Create(
+      sparse_linear_algebra_library_type, ordering_type));
+  const CompressedRowSparseMatrix::StorageType storage_type =
+      sparse_cholesky->StorageType();
+
+  scoped_ptr<CompressedRowSparseMatrix> lhs(
+      CreateRandomSymmetricPositiveDefiniteMatrix(num_blocks,
+                                                  min_block_size,
+                                                  max_block_size,
+                                                  block_density,
+                                                  storage_type));
+  if (!use_block_structure) {
+    lhs->mutable_row_blocks()->clear();
+    lhs->mutable_col_blocks()->clear();
+  }
+
+  Vector rhs = Vector::Random(lhs->num_rows());
+  Vector expected(lhs->num_rows());
+  Vector actual(lhs->num_rows());
+
+  EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));
+  std::string message;
+  EXPECT_EQ(sparse_cholesky->FactorAndSolve(
+                lhs.get(), rhs.data(), actual.data(), &message),
+            LINEAR_SOLVER_SUCCESS);
+  Matrix eigen_lhs;
+  lhs->ToDenseMatrix(&eigen_lhs);
+  EXPECT_NEAR((actual - expected).norm() / actual.norm(),
+              0.0,
+              std::numeric_limits<double>::epsilon() * 10)
+      << "\n"
+      << eigen_lhs;
+}
+
+typedef ::std::tr1::tuple<SparseLinearAlgebraLibraryType, OrderingType, bool>
+    Param;
+
+std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+  Param param = info.param;
+  std::stringstream ss;
+  ss << SparseLinearAlgebraLibraryTypeToString(std::tr1::get<0>(param)) << "_"
+     << (std::tr1::get<1>(param) == AMD ? "AMD" : "NATURAL") << "_"
+     << (std::tr1::get<2>(param) ? "UseBlockStructure" : "NoBlockStructure");
+  return ss.str();
+}
+
+class SparseCholeskyTest : public ::testing::TestWithParam<Param> {};
+
+TEST_P(SparseCholeskyTest, FactorAndSolve) {
+  SetRandomState(2982);
+  const int kMinNumBlocks = 1;
+  const int kMaxNumBlocks = 10;
+  const int kNumTrials = 10;
+  const int kMinBlockSize = 1;
+  const int kMaxBlockSize = 5;
+
+  for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
+       ++num_blocks) {
+    for (int trial = 0; trial < kNumTrials; ++trial) {
+      const double block_density = std::max(0.1, RandDouble());
+      Param param = GetParam();
+      SparseCholeskySolverUnitTest(std::tr1::get<0>(param),
+                                   std::tr1::get<1>(param),
+                                   std::tr1::get<2>(param),
+                                   num_blocks,
+                                   kMinBlockSize,
+                                   kMaxBlockSize,
+                                   block_density);
+    }
+  }
+}
+
+#ifndef CERES_NO_SUITESPARSE
+INSTANTIATE_TEST_CASE_P(SuiteSparseCholesky,
+                        SparseCholeskyTest,
+                        ::testing::Combine(::testing::Values(SUITE_SPARSE),
+                                           ::testing::Values(AMD, NATURAL),
+                                           ::testing::Values(true, false)),
+                        ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_CXSPARSE
+INSTANTIATE_TEST_CASE_P(CXSparseCholesky,
+                        SparseCholeskyTest,
+                        ::testing::Combine(::testing::Values(CX_SPARSE),
+                                           ::testing::Values(AMD, NATURAL),
+                                           ::testing::Values(true, false)),
+                        ParamInfoToString);
+#endif
+
+#ifdef CERES_USE_EIGEN_SPARSE
+INSTANTIATE_TEST_CASE_P(EigenSparseCholesky,
+                        SparseCholeskyTest,
+                        ::testing::Combine(::testing::Values(EIGEN_SPARSE),
+                                           ::testing::Values(AMD, NATURAL),
+                                           ::testing::Values(true, false)),
+                        ParamInfoToString);
+#endif
+
+}  // namespace internal
+}  // namespace ceres

internal/ceres/sparse_normal_cholesky_solver.cc

@@ -33,138 +33,44 @@
 #include <algorithm>
 #include <cstring>
 #include <ctime>
-#include <sstream>
 
-#include "Eigen/SparseCore"
 #include "ceres/compressed_row_sparse_matrix.h"
-#include "ceres/cxsparse.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/linear_solver.h"
-#include "ceres/suitesparse.h"
+#include "ceres/sparse_cholesky.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "ceres/types.h"
 #include "ceres/wall_time.h"
 
-#ifdef CERES_USE_EIGEN_SPARSE
-#include "Eigen/SparseCholesky"
-#endif
-
 namespace ceres {
 namespace internal {
 
-// Different sparse linear algebra libraries prefer different storage
-// orders for the input matrix. This trait class helps choose the
-// ordering based on the sparse linear algebra backend being used.
-//
-// The storage order is lower-triangular by default. It is only
-// SuiteSparse which prefers an upper triangular matrix. Saves a whole
-// matrix copy in the process.
-//
-// Note that this is the storage order for a compressed row sparse
-// matrix. All the sparse linear algebra libraries take compressed
-// column sparse matrices as input. We map these matrices to into
-// compressed column sparse matrices before calling them and in the
-// process, transpose them.
-//
-// TODO(sameeragarwal): This does not account for post ordering, where
-// the optimal storage order maybe different. Either get rid of post
-// ordering support entirely, or investigate making this trait class
-// richer.
-
-CompressedRowSparseMatrix::StorageType StorageTypeForSparseLinearAlgebraLibrary(
-    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type) {
-  if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
-    return CompressedRowSparseMatrix::UPPER_TRIANGULAR;
-  }
-  return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
-}
-
-namespace {
-
-#ifdef CERES_USE_EIGEN_SPARSE
-// A templated factorized and solve function, which allows us to use
-// the same code independent of whether a AMD or a Natural ordering is
-// used.
-template <typename SimplicialCholeskySolver, typename SparseMatrixType>
-LinearSolver::Summary SimplicialLDLTSolve(const SparseMatrixType& lhs,
-                                          const bool do_symbolic_analysis,
-                                          SimplicialCholeskySolver* solver,
-                                          double* rhs_and_solution,
-                                          EventLogger* event_logger) {
-  LinearSolver::Summary summary;
-  summary.num_iterations = 1;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  summary.message = "Success.";
-
-  if (do_symbolic_analysis) {
-    solver->analyzePattern(lhs);
-    if (VLOG_IS_ON(2)) {
-      std::stringstream ss;
-      solver->dumpMemory(ss);
-      VLOG(2) << "Symbolic Analysis\n" << ss.str();
-    }
-    event_logger->AddEvent("Analyze");
-    if (solver->info() != Eigen::Success) {
-      summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-      summary.message = "Eigen failure. Unable to find symbolic factorization.";
-      return summary;
-    }
-  }
-
-  solver->factorize(lhs);
-  event_logger->AddEvent("Factorize");
-  if (solver->info() != Eigen::Success) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "Eigen failure. Unable to find numeric factorization.";
-    return summary;
-  }
-
-  const Vector rhs = VectorRef(rhs_and_solution, lhs.cols());
-
-  VectorRef(rhs_and_solution, lhs.cols()) = solver->solve(rhs);
-  event_logger->AddEvent("Solve");
-  if (solver->info() != Eigen::Success) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "Eigen failure. Unable to do triangular solve.";
-    return summary;
-  }
-
-  return summary;
-}
-
-#endif  // CERES_USE_EIGEN_SPARSE
-
-}  // namespace
-
 SparseNormalCholeskySolver::SparseNormalCholeskySolver(
     const LinearSolver::Options& options)
-    : factor_(NULL), cxsparse_factor_(NULL), options_(options) {}
-
-void SparseNormalCholeskySolver::FreeFactorization() {
-  if (factor_ != NULL) {
-    ss_.Free(factor_);
-    factor_ = NULL;
-  }
-
-  if (cxsparse_factor_ != NULL) {
-    cxsparse_.Free(cxsparse_factor_);
-    cxsparse_factor_ = NULL;
-  }
+    : options_(options) {
+  sparse_cholesky_.reset(
+      SparseCholesky::Create(options_.sparse_linear_algebra_library_type,
+                             options_.use_postordering ? AMD : NATURAL));
 }
 
-SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {
-  FreeFactorization();
-}
+SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {}
 
 LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
     CompressedRowSparseMatrix* A,
     const double* b,
     const LinearSolver::PerSolveOptions& per_solve_options,
     double* x) {
+  EventLogger event_logger("SparseNormalCholeskySolver::Solve");
+  LinearSolver::Summary summary;
+  summary.num_iterations = 1;
+  summary.termination_type = LINEAR_SOLVER_SUCCESS;
+  summary.message = "Success.";
+
   const int num_cols = A->num_cols();
   VectorRef(x, num_cols).setZero();
   A->LeftMultiply(b, x);
+  event_logger.AddEvent("Compute RHS");
 
   if (per_solve_options.D != NULL) {
     // Temporarily append a diagonal block to the A matrix, but undo
@@ -179,241 +85,28 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
     }
     A->AppendRows(*regularizer);
   }
+  event_logger.AddEvent("Append Rows");
 
   if (outer_product_.get() == NULL) {
     outer_product_.reset(
         CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
-            *A,
-            StorageTypeForSparseLinearAlgebraLibrary(
-                options_.sparse_linear_algebra_library_type),
-            &pattern_));
+            *A, sparse_cholesky_->StorageType(), &pattern_));
+    event_logger.AddEvent("Outer Product Program");
   }
 
   CompressedRowSparseMatrix::ComputeOuterProduct(
       *A, pattern_, outer_product_.get());
-
-  LinearSolver::Summary summary;
-  switch (options_.sparse_linear_algebra_library_type) {
-    case SUITE_SPARSE:
-      summary = SolveImplUsingSuiteSparse(x);
-      break;
-    case CX_SPARSE:
-      summary = SolveImplUsingCXSparse(x);
-      break;
-    case EIGEN_SPARSE:
-      summary = SolveImplUsingEigen(x);
-      break;
-    default:
-      LOG(FATAL) << "Unknown sparse linear algebra library : "
-                 << options_.sparse_linear_algebra_library_type;
-  }
+  event_logger.AddEvent("Outer Product");
 
   if (per_solve_options.D != NULL) {
     A->DeleteRows(num_cols);
   }
 
+  summary.termination_type = sparse_cholesky_->FactorAndSolve(
+      outer_product_.get(), x, x, &summary.message);
+  event_logger.AddEvent("Factor & Solve");
   return summary;
 }
 
-LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingEigen(
-    double* rhs_and_solution) {
-#ifndef CERES_USE_EIGEN_SPARSE
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-  summary.message =
-      "SPARSE_NORMAL_CHOLESKY cannot be used with EIGEN_SPARSE "
-      "because Ceres was not built with support for "
-      "Eigen's SimplicialLDLT decomposition. "
-      "This requires enabling building with -DEIGENSPARSE=ON.";
-  return summary;
-
-#else
-
-  EventLogger event_logger("SparseNormalCholeskySolver::Eigen::Solve");
-
-  // Map outer_product_ to an upper triangular column major matrix.
-  //
-  // 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.
-  Eigen::MappedSparseMatrix<double, Eigen::ColMajor> lhs(
-      outer_product_->num_rows(),
-      outer_product_->num_rows(),
-      outer_product_->num_nonzeros(),
-      outer_product_->mutable_rows(),
-      outer_product_->mutable_cols(),
-      outer_product_->mutable_values());
-
-  bool do_symbolic_analysis = false;
-
-  // If using post ordering or an old version of Eigen, we cannot
-  // depend on a preordered jacobian, so we work with a SimplicialLDLT
-  // decomposition with AMD ordering.
-  if (options_.use_postordering || !EIGEN_VERSION_AT_LEAST(3, 2, 2)) {
-    if (amd_ldlt_.get() == NULL) {
-      amd_ldlt_.reset(new SimplicialLDLTWithAMDOrdering);
-      do_symbolic_analysis = true;
-    }
-
-    return SimplicialLDLTSolve(lhs,
-                               do_symbolic_analysis,
-                               amd_ldlt_.get(),
-                               rhs_and_solution,
-                               &event_logger);
-  }
-
-#if EIGEN_VERSION_AT_LEAST(3, 2, 2)
-  // The common case
-  if (natural_ldlt_.get() == NULL) {
-    natural_ldlt_.reset(new SimplicialLDLTWithNaturalOrdering);
-    do_symbolic_analysis = true;
-  }
-
-  return SimplicialLDLTSolve(lhs,
-                             do_symbolic_analysis,
-                             natural_ldlt_.get(),
-                             rhs_and_solution,
-                             &event_logger);
-#endif
-
-#endif  // EIGEN_USE_EIGEN_SPARSE
-}
-
-LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
-    double* rhs_and_solution) {
-#ifdef CERES_NO_CXSPARSE
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-  summary.message =
-      "SPARSE_NORMAL_CHOLESKY cannot be used with CX_SPARSE "
-      "because Ceres was not built with support for CXSparse. "
-      "This requires enabling building with -DCXSPARSE=ON.";
-
-  return summary;
-
-#else
-
-  EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve");
-  LinearSolver::Summary summary;
-  summary.num_iterations = 1;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  summary.message = "Success.";
-
-  // Map outer_product_ to an upper triangular column major matrix.
-  //
-  // 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.
-  cs_di lhs = cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
-
-  event_logger.AddEvent("Setup");
-
-  // Compute symbolic factorization if not available.
-  if (cxsparse_factor_ == NULL) {
-    if (options_.use_postordering) {
-      cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(
-          &lhs, outer_product_->col_blocks(), outer_product_->col_blocks());
-    } else {
-      cxsparse_factor_ = cxsparse_.AnalyzeCholeskyWithNaturalOrdering(&lhs);
-    }
-  }
-  event_logger.AddEvent("Analysis");
-
-  if (cxsparse_factor_ == NULL) {
-    summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-    summary.message =
-        "CXSparse failure. Unable to find symbolic factorization.";
-  } else if (!cxsparse_.SolveCholesky(
-                 &lhs, cxsparse_factor_, rhs_and_solution)) {
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-    summary.message = "CXSparse::SolveCholesky failed.";
-  }
-  event_logger.AddEvent("Solve");
-
-  return summary;
-#endif
-}
-
-LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
-    double* rhs_and_solution) {
-#ifdef CERES_NO_SUITESPARSE
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 0;
-  summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-  summary.message =
-      "SPARSE_NORMAL_CHOLESKY cannot be used with SUITE_SPARSE "
-      "because Ceres was not built with support for SuiteSparse. "
-      "This requires enabling building with -DSUITESPARSE=ON.";
-  return summary;
-
-#else
-
-  EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve");
-  LinearSolver::Summary summary;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  summary.num_iterations = 1;
-  summary.message = "Success.";
-
-  // Map outer_product_ to a lower triangular column major matrix.
-  //
-  // outer_product_ is a compressed row sparse matrix and in upper
-  // triangular form, when mapped to a compressed column sparse
-  // matrix, it becomes a lower triangular matrix.
-  const int num_cols = outer_product_->num_cols();
-  cholmod_sparse lhs =
-      ss_.CreateSparseMatrixTransposeView(outer_product_.get());
-  event_logger.AddEvent("Setup");
-
-  if (factor_ == NULL) {
-    if (options_.use_postordering) {
-      factor_ = ss_.BlockAnalyzeCholesky(
-          &lhs,
-          outer_product_->col_blocks(),
-          outer_product_->col_blocks(),
-          &summary.message);
-    } else {
-      factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message);
-    }
-  }
-
-  event_logger.AddEvent("Analysis");
-
-  if (factor_ == NULL) {
-    summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
-    // No need to set message as it has already been set by the
-    // symbolic analysis routines above.
-    return summary;
-  }
-
-  summary.termination_type = ss_.Cholesky(&lhs, factor_, &summary.message);
-  if (summary.termination_type != LINEAR_SOLVER_SUCCESS) {
-    return summary;
-  }
-
-  cholmod_dense* rhs =
-      ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols);
-  cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message);
-  event_logger.AddEvent("Solve");
-
-  ss_.Free(rhs);
-  if (solution != NULL) {
-    memcpy(rhs_and_solution, solution->x, num_cols * sizeof(*rhs_and_solution));
-    ss_.Free(solution);
-  } else {
-    // No need to set message as it has already been set by the
-    // numeric factorization routine above.
-    summary.termination_type = LINEAR_SOLVER_FAILURE;
-  }
-
-  event_logger.AddEvent("Teardown");
-  return summary;
-#endif
-}
-
 }  // namespace internal
 }  // namespace ceres

internal/ceres/sparse_normal_cholesky_solver.h

@@ -34,27 +34,21 @@
 #ifndef CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
 #define CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
 
-#include <vector>
-
 // This include must come before any #ifndef check on Ceres compile options.
 #include "ceres/internal/port.h"
 
+#include <vector>
 #include "ceres/internal/macros.h"
 #include "ceres/linear_solver.h"
-#include "ceres/suitesparse.h"
-#include "ceres/cxsparse.h"
-
-#ifdef CERES_USE_EIGEN_SPARSE
-#include "Eigen/SparseCholesky"
-#endif
 
 namespace ceres {
 namespace internal {
 
 class CompressedRowSparseMatrix;
+class SparseCholesky;
 
-// Solves the normal equations (A'A + D'D) x = A'b, using the CHOLMOD sparse
-// cholesky solver.
+// Solves the normal equations (A'A + D'D) x = A'b, using the sparse
+// linear algebra library of the user's choice.
 class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver {
  public:
   explicit SparseNormalCholeskySolver(const LinearSolver::Options& options);
@@ -67,48 +61,10 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver {
       const LinearSolver::PerSolveOptions& options,
       double* x);
 
-  LinearSolver::Summary SolveImplUsingSuiteSparse(double* rhs_and_solution);
-  LinearSolver::Summary SolveImplUsingCXSparse(double* rhs_and_solution);
-  LinearSolver::Summary SolveImplUsingEigen(double* rhs_and_solution);
-
-  void FreeFactorization();
-
-  SuiteSparse ss_;
-  // Cached factorization
-  cholmod_factor* factor_;
-
-  CXSparse cxsparse_;
-  // Cached factorization
-  cs_dis* cxsparse_factor_;
-
-#ifdef CERES_USE_EIGEN_SPARSE
-
-  // The preprocessor gymnastics here are dealing with the fact that
-  // before version 3.2.2, Eigen did not support a third template
-  // parameter to specify the ordering.
-#if EIGEN_VERSION_AT_LEAST(3,2,2)
-  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper,
-                                Eigen::NaturalOrdering<int> >
-  SimplicialLDLTWithNaturalOrdering;
-  scoped_ptr<SimplicialLDLTWithNaturalOrdering> natural_ldlt_;
-
-  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper,
-                                Eigen::AMDOrdering<int> >
-  SimplicialLDLTWithAMDOrdering;
-  scoped_ptr<SimplicialLDLTWithAMDOrdering> amd_ldlt_;
-
-#else
-  typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper>
-  SimplicialLDLTWithAMDOrdering;
-
-  scoped_ptr<SimplicialLDLTWithAMDOrdering> amd_ldlt_;
-#endif
-
-#endif
-
+  const LinearSolver::Options options_;
+  scoped_ptr<SparseCholesky> sparse_cholesky_;
   scoped_ptr<CompressedRowSparseMatrix> outer_product_;
   std::vector<int> pattern_;
-  const LinearSolver::Options options_;
   CERES_DISALLOW_COPY_AND_ASSIGN(SparseNormalCholeskySolver);
 };
 

internal/ceres/suitesparse.cc

@@ -35,11 +35,12 @@
 #include "ceres/suitesparse.h"
 
 #include <vector>
-#include "cholmod.h"
+
 #include "ceres/compressed_col_sparse_matrix_utils.h"
 #include "ceres/compressed_row_sparse_matrix.h"
 #include "ceres/linear_solver.h"
 #include "ceres/triplet_sparse_matrix.h"
+#include "cholmod.h"
 
 namespace ceres {
 namespace internal {
@@ -47,13 +48,9 @@ namespace internal {
 using std::string;
 using std::vector;
 
-SuiteSparse::SuiteSparse() {
-  cholmod_start(&cc_);
-}
+SuiteSparse::SuiteSparse() { cholmod_start(&cc_); }
 
-SuiteSparse::~SuiteSparse() {
-  cholmod_finish(&cc_);
-}
+SuiteSparse::~SuiteSparse() { cholmod_finish(&cc_); }
 
 cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) {
   cholmod_triplet triplet;
@@ -73,7 +70,6 @@ cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) {
   return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_);
 }
 
-
 cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose(
     TripletSparseMatrix* A) {
   cholmod_triplet triplet;
@@ -127,12 +123,12 @@ cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView(
 cholmod_dense* SuiteSparse::CreateDenseVector(const double* x,
                                               int in_size,
                                               int out_size) {
-    CHECK_LE(in_size, out_size);
-    cholmod_dense* v = cholmod_zeros(out_size, 1, CHOLMOD_REAL, &cc_);
-    if (x != NULL) {
-      memcpy(v->x, x, in_size*sizeof(*x));
-    }
-    return v;
+  CHECK_LE(in_size, out_size);
+  cholmod_dense* v = cholmod_zeros(out_size, 1, CHOLMOD_REAL, &cc_);
+  if (x != NULL) {
+    memcpy(v->x, x, in_size * sizeof(*x));
+  }
+  return v;
 }
 
 cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A,
@@ -151,19 +147,18 @@ cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A,
   }
 
   if (cc_.status != CHOLMOD_OK) {
-    *message = StringPrintf("cholmod_analyze failed. error code: %d",
-                            cc_.status);
+    *message =
+        StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
     return NULL;
   }
 
   return CHECK_NOTNULL(factor);
 }
 
-cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(
-    cholmod_sparse* A,
-    const vector<int>& row_blocks,
-    const vector<int>& col_blocks,
-    string* message) {
+cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(cholmod_sparse* A,
+                                                  const vector<int>& row_blocks,
+                                                  const vector<int>& col_blocks,
+                                                  string* message) {
   vector<int> ordering;
   if (!BlockAMDOrdering(A, row_blocks, col_blocks, &ordering)) {
     return NULL;
@@ -172,22 +167,20 @@ cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(
 }
 
 cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(
-    cholmod_sparse* A,
-    const vector<int>& ordering,
-    string* message) {
+    cholmod_sparse* A, const vector<int>& ordering, string* message) {
   CHECK_EQ(ordering.size(), A->nrow);
 
   cc_.nmethods = 1;
   cc_.method[0].ordering = CHOLMOD_GIVEN;
 
-  cholmod_factor* factor  =
+  cholmod_factor* factor =
       cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), NULL, 0, &cc_);
   if (VLOG_IS_ON(2)) {
     cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
   }
   if (cc_.status != CHOLMOD_OK) {
-    *message = StringPrintf("cholmod_analyze failed. error code: %d",
-                            cc_.status);
+    *message =
+        StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
     return NULL;
   }
 
@@ -195,19 +188,18 @@ cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(
 }
 
 cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering(
-    cholmod_sparse* A,
-    string* message) {
+    cholmod_sparse* A, string* message) {
   cc_.nmethods = 1;
   cc_.method[0].ordering = CHOLMOD_NATURAL;
   cc_.postorder = 0;
 
-  cholmod_factor* factor  = cholmod_analyze(A, &cc_);
+  cholmod_factor* factor = cholmod_analyze(A, &cc_);
   if (VLOG_IS_ON(2)) {
     cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
   }
   if (cc_.status != CHOLMOD_OK) {
-    *message = StringPrintf("cholmod_analyze failed. error code: %d",
-                            cc_.status);
+    *message =
+        StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
     return NULL;
   }
 
@@ -232,7 +224,6 @@ bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A,
                                             col_blocks,
                                             &block_rows,
                                             &block_cols);
-
   cholmod_sparse_struct block_matrix;
   block_matrix.nrow = num_row_blocks;
   block_matrix.ncol = num_col_blocks;
@@ -273,13 +264,6 @@ LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A,
   int cholmod_status = cholmod_factorize(A, L, &cc_);
   cc_.print = old_print_level;
 
-  // TODO(sameeragarwal): This switch statement is not consistent. It
-  // treats all kinds of CHOLMOD failures as warnings. Some of these
-  // like out of memory are definitely not warnings. The problem is
-  // that the return value Cholesky is two valued, but the state of
-  // the linear solver is really three valued. SUCCESS,
-  // NON_FATAL_FAILURE (e.g., indefinite matrix) and FATAL_FAILURE
-  // (e.g. out of memory).
   switch (cc_.status) {
     case CHOLMOD_NOT_INSTALLED:
       *message = "CHOLMOD failure: Method not installed.";
@@ -297,23 +281,25 @@ LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A,
       *message = "CHOLMOD warning: Matrix not positive definite.";
       return LINEAR_SOLVER_FAILURE;
     case CHOLMOD_DSMALL:
-      *message = "CHOLMOD warning: D for LDL' or diag(L) or "
-                "LL' has tiny absolute value.";
+      *message =
+          "CHOLMOD warning: D for LDL' or diag(L) or "
+          "LL' has tiny absolute value.";
       return LINEAR_SOLVER_FAILURE;
     case CHOLMOD_OK:
       if (cholmod_status != 0) {
         return LINEAR_SOLVER_SUCCESS;
       }
 
-      *message = "CHOLMOD failure: cholmod_factorize returned false "
+      *message =
+          "CHOLMOD failure: cholmod_factorize returned false "
           "but cholmod_common::status is CHOLMOD_OK."
           "Please report this to ceres-solver@googlegroups.com.";
       return LINEAR_SOLVER_FATAL_ERROR;
     default:
-      *message =
-          StringPrintf("Unknown cholmod return code: %d. "
-                       "Please report this to ceres-solver@googlegroups.com.",
-                       cc_.status);
+      *message = StringPrintf(
+          "Unknown cholmod return code: %d. "
+          "Please report this to ceres-solver@googlegroups.com.",
+          cc_.status);
       return LINEAR_SOLVER_FATAL_ERROR;
   }
 
@@ -337,9 +323,7 @@ bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
 }
 
 bool SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering(
-    cholmod_sparse* matrix,
-    int* constraints,
-    int* ordering) {
+    cholmod_sparse* matrix, int* constraints, int* ordering) {
 #ifndef CERES_NO_CAMD
   return cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_);
 #else
@@ -352,6 +336,79 @@ bool SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering(
 #endif
 }
 
+SuiteSparseCholesky* SuiteSparseCholesky::Create(
+    const OrderingType ordering_type) {
+  return new SuiteSparseCholesky(ordering_type);
+}
+
+SuiteSparseCholesky::SuiteSparseCholesky(const OrderingType ordering_type)
+    : ordering_type_(ordering_type), factor_(NULL) {}
+
+SuiteSparseCholesky::~SuiteSparseCholesky() {
+  if (factor_ != NULL) {
+    ss_.Free(factor_);
+  }
+}
+
+LinearSolverTerminationType SuiteSparseCholesky::Factorize(
+    CompressedRowSparseMatrix* lhs, string* message) {
+  if (lhs == NULL) {
+    *message = "Failure: Input lhs is NULL.";
+    return LINEAR_SOLVER_FATAL_ERROR;
+  }
+
+  cholmod_sparse cholmod_lhs = ss_.CreateSparseMatrixTransposeView(lhs);
+
+  if (factor_ == NULL) {
+    if (ordering_type_ == NATURAL) {
+      factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&cholmod_lhs, message);
+    } else {
+      if (!lhs->col_blocks().empty() && !(lhs->row_blocks().empty())) {
+        factor_ = ss_.BlockAnalyzeCholesky(
+            &cholmod_lhs, lhs->col_blocks(), lhs->row_blocks(), message);
+      } else {
+        factor_ = ss_.AnalyzeCholesky(&cholmod_lhs, message);
+      }
+    }
+
+    if (factor_ == NULL) {
+      return LINEAR_SOLVER_FATAL_ERROR;
+    }
+  }
+
+  return ss_.Cholesky(&cholmod_lhs, factor_, message);
+}
+
+CompressedRowSparseMatrix::StorageType SuiteSparseCholesky::StorageType()
+    const {
+  return ((ordering_type_ == NATURAL)
+              ? CompressedRowSparseMatrix::UPPER_TRIANGULAR
+              : CompressedRowSparseMatrix::LOWER_TRIANGULAR);
+}
+
+LinearSolverTerminationType SuiteSparseCholesky::Solve(const double* rhs,
+                                                       double* solution,
+                                                       string* message) {
+  // Error checking
+  if (factor_ == NULL) {
+    *message = "Solve called without a call to Factorize first.";
+    return LINEAR_SOLVER_FATAL_ERROR;
+  }
+
+  const int num_cols = factor_->n;
+  cholmod_dense* cholmod_dense_rhs =
+      ss_.CreateDenseVector(rhs, num_cols, num_cols);
+  cholmod_dense* cholmod_dense_solution =
+      ss_.Solve(factor_, cholmod_dense_rhs, message);
+  ss_.Free(cholmod_dense_rhs);
+  if (cholmod_dense_solution == NULL) {
+    return LINEAR_SOLVER_FAILURE;
+  }
+
+  memcpy(solution, cholmod_dense_solution->x, num_cols * sizeof(*solution));
+  return LINEAR_SOLVER_SUCCESS;
+}
+
 }  // namespace internal
 }  // namespace ceres
 

internal/ceres/suitesparse.h

@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2017 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -41,11 +41,11 @@
 #include <cstring>
 #include <string>
 #include <vector>
-
+#include "SuiteSparseQR.hpp"
 #include "ceres/linear_solver.h"
+#include "ceres/sparse_cholesky.h"
 #include "cholmod.h"
 #include "glog/logging.h"
-#include "SuiteSparseQR.hpp"
 
 // Before SuiteSparse version 4.2.0, cholmod_camd was only enabled
 // if SuiteSparse was compiled with Metis support. This makes
@@ -278,6 +278,27 @@ class SuiteSparse {
   cholmod_common cc_;
 };
 
+class SuiteSparseCholesky : public SparseCholesky {
+ public:
+  static SuiteSparseCholesky* Create(const OrderingType ordering_type);
+
+  // SparseCholesky interface.
+  virtual ~SuiteSparseCholesky();
+  virtual CompressedRowSparseMatrix::StorageType StorageType() const;
+  virtual LinearSolverTerminationType Factorize(
+      CompressedRowSparseMatrix* lhs, std::string* message);
+  virtual LinearSolverTerminationType Solve(const double* rhs,
+                                            double* solution,
+                                            std::string* message);
+ private:
+  SuiteSparseCholesky(const OrderingType ordering_type);
+
+  const OrderingType ordering_type_;
+  SuiteSparse ss_;
+  cholmod_factor* factor_;
+  cholmod_dense* rhs_;
+};
+
 }  // namespace internal
 }  // namespace ceres
 
@@ -285,6 +306,9 @@ class SuiteSparse {
 
 typedef void cholmod_factor;
 
+namespace ceres {
+namespace internal {
+
 class SuiteSparse {
  public:
   // Defining this static function even when SuiteSparse is not
@@ -301,6 +325,9 @@ class SuiteSparse {
   void Free(void* arg) {}
 };
 
+}  // namespace internal
+}  // namespace ceres
+
 #endif  // CERES_NO_SUITESPARSE
 
 #endif  // CERES_INTERNAL_SUITESPARSE_H_

internal/ceres/unsymmetric_linear_solver_test.cc

@@ -38,10 +38,10 @@
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
-
 namespace ceres {
 namespace internal {
 
+// TODO(sameeragarwal): Refactor and expand these tests.
 class UnsymmetricLinearSolverTest : public ::testing::Test {
  protected :
   virtual void SetUp() {
@@ -254,7 +254,6 @@ TEST_F(UnsymmetricLinearSolverTest,
   options.dynamic_sparsity = true;
   TestSolver(options);
 }
-
 #endif  // CERES_USE_EIGEN_SPARSE
 
 }  // namespace internal

internal/ceres/visibility.cc

@@ -28,11 +28,6 @@
 //
 // Author: kushalav@google.com (Avanish Kushal)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/visibility.h"
 
 #include <cmath>
@@ -162,5 +157,3 @@ WeightedGraph<int>* CreateSchurComplementGraph(
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/visibility.h

@@ -35,11 +35,6 @@
 #ifndef CERES_INTERNAL_VISIBILITY_H_
 #define CERES_INTERNAL_VISIBILITY_H_
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include <set>
 #include <vector>
 #include "ceres/graph.h"
@@ -80,5 +75,4 @@ WeightedGraph<int>* CreateSchurComplementGraph(
 }  // namespace internal
 }  // namespace ceres
 
-#endif  // CERES_NO_SUITESPARSE
 #endif  // CERES_INTERNAL_VISIBILITY_H_

internal/ceres/visibility_based_preconditioner.cc

@@ -28,11 +28,6 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/visibility_based_preconditioner.h"
 
 #include <algorithm>
@@ -77,19 +72,14 @@ static const double kSingleLinkageMinSimilarity = 0.9;
 VisibilityBasedPreconditioner::VisibilityBasedPreconditioner(
     const CompressedRowBlockStructure& bs,
     const Preconditioner::Options& options)
-    : options_(options),
-      num_blocks_(0),
-      num_clusters_(0),
-      factor_(NULL) {
+    : options_(options), num_blocks_(0), num_clusters_(0) {
   CHECK_GT(options_.elimination_groups.size(), 1);
   CHECK_GT(options_.elimination_groups[0], 0);
-  CHECK(options_.type == CLUSTER_JACOBI ||
-        options_.type == CLUSTER_TRIDIAGONAL)
+  CHECK(options_.type == CLUSTER_JACOBI || options_.type == CLUSTER_TRIDIAGONAL)
       << "Unknown preconditioner type: " << options_.type;
   num_blocks_ = bs.cols.size() - options_.elimination_groups[0];
-  CHECK_GT(num_blocks_, 0)
-      << "Jacobian should have atleast 1 f_block for "
-      << "visibility based preconditioning.";
+  CHECK_GT(num_blocks_, 0) << "Jacobian should have atleast 1 f_block for "
+                           << "visibility based preconditioning.";
 
   // Vector of camera block sizes
   block_size_.resize(num_blocks_);
@@ -114,29 +104,17 @@ VisibilityBasedPreconditioner::VisibilityBasedPreconditioner(
   InitEliminator(bs);
   const time_t eliminator_time = time(NULL);
 
-  // Allocate temporary storage for a vector used during
-  // RightMultiply.
-  tmp_rhs_ = CHECK_NOTNULL(ss_.CreateDenseVector(NULL,
-                                                 m_->num_rows(),
-                                                 m_->num_rows()));
+  sparse_cholesky_.reset(
+      SparseCholesky::Create(options_.sparse_linear_algebra_library_type, AMD));
+
   const time_t init_time = time(NULL);
-  VLOG(2) << "init time: "
-          << init_time - start_time
+  VLOG(2) << "init time: " << init_time - start_time
           << " structure time: " << structure_time - start_time
           << " storage time:" << storage_time - structure_time
           << " eliminator time: " << eliminator_time - storage_time;
 }
 
-VisibilityBasedPreconditioner::~VisibilityBasedPreconditioner() {
-  if (factor_ != NULL) {
-    ss_.Free(factor_);
-    factor_ = NULL;
-  }
-  if (tmp_rhs_ != NULL) {
-    ss_.Free(tmp_rhs_);
-    tmp_rhs_ = NULL;
-  }
-}
+VisibilityBasedPreconditioner::~VisibilityBasedPreconditioner() {}
 
 // Determine the sparsity structure of the CLUSTER_JACOBI
 // preconditioner. It clusters cameras using their scene
@@ -203,22 +181,17 @@ void VisibilityBasedPreconditioner::ClusterCameras(
   if (options_.visibility_clustering_type == CANONICAL_VIEWS) {
     vector<int> centers;
     CanonicalViewsClusteringOptions clustering_options;
-    clustering_options.size_penalty_weight =
-        kCanonicalViewsSizePenaltyWeight;
+    clustering_options.size_penalty_weight = kCanonicalViewsSizePenaltyWeight;
     clustering_options.similarity_penalty_weight =
         kCanonicalViewsSimilarityPenaltyWeight;
-    ComputeCanonicalViewsClustering(clustering_options,
-                                    *schur_complement_graph,
-                                    &centers,
-                                    &membership);
+    ComputeCanonicalViewsClustering(
+        clustering_options, *schur_complement_graph, &centers, &membership);
     num_clusters_ = centers.size();
   } else if (options_.visibility_clustering_type == SINGLE_LINKAGE) {
     SingleLinkageClusteringOptions clustering_options;
-    clustering_options.min_similarity =
-        kSingleLinkageMinSimilarity;
-    num_clusters_ = ComputeSingleLinkageClustering(clustering_options,
-                                                   *schur_complement_graph,
-                                                   &membership);
+    clustering_options.min_similarity = kSingleLinkageMinSimilarity;
+    num_clusters_ = ComputeSingleLinkageClustering(
+        clustering_options, *schur_complement_graph, &membership);
   } else {
     LOG(FATAL) << "Unknown visibility clustering algorithm.";
   }
@@ -417,13 +390,12 @@ void VisibilityBasedPreconditioner::ScaleOffDiagonalCells() {
     }
 
     int r, c, row_stride, col_stride;
-    CellInfo* cell_info = m_->GetCell(block1, block2,
-                                      &r, &c,
-                                      &row_stride, &col_stride);
+    CellInfo* cell_info =
+        m_->GetCell(block1, block2, &r, &c, &row_stride, &col_stride);
     CHECK(cell_info != NULL)
         << "Cell missing for block pair (" << block1 << "," << block2 << ")"
-        << " cluster pair (" << cluster_membership_[block1]
-        << " " << cluster_membership_[block2] << ")";
+        << " cluster pair (" << cluster_membership_[block1] << " "
+        << cluster_membership_[block2] << ")";
 
     // Ah the magic of tri-diagonal matrices and diagonal
     // dominance. See Lemma 1 in "Visibility Based Preconditioning
@@ -437,58 +409,42 @@ void VisibilityBasedPreconditioner::ScaleOffDiagonalCells() {
 // matrix.
 LinearSolverTerminationType VisibilityBasedPreconditioner::Factorize() {
   // Extract the TripletSparseMatrix that is used for actually storing
-  // S and convert it into a cholmod_sparse object.
-  cholmod_sparse* lhs = ss_.CreateSparseMatrix(
-      down_cast<BlockRandomAccessSparseMatrix*>(
-          m_.get())->mutable_matrix());
-
-  // The matrix is symmetric, and the upper triangular part of the
-  // matrix contains the values.
-  lhs->stype = 1;
-
-  // TODO(sameeragarwal): Refactor to pipe this up and out.
-  std::string status;
-
-  // Symbolic factorization is computed if we don't already have one handy.
-  if (factor_ == NULL) {
-    factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_, &status);
+  // S and convert it into a CompressedRowSparseMatrix.
+  const TripletSparseMatrix* tsm =
+      down_cast<BlockRandomAccessSparseMatrix*>(m_.get())->mutable_matrix();
+
+  scoped_ptr<CompressedRowSparseMatrix> lhs;
+  const CompressedRowSparseMatrix::StorageType storage_type =
+      sparse_cholesky_->StorageType();
+  if (storage_type == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+    lhs.reset(CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
+    lhs->set_storage_type(CompressedRowSparseMatrix::UPPER_TRIANGULAR);
+  } else {
+    lhs.reset(
+        CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm));
+    lhs->set_storage_type(CompressedRowSparseMatrix::LOWER_TRIANGULAR);
   }
 
-  const LinearSolverTerminationType termination_type =
-      (factor_ != NULL)
-      ? ss_.Cholesky(lhs, factor_, &status)
-      : LINEAR_SOLVER_FATAL_ERROR;
-
-  ss_.Free(lhs);
-  return termination_type;
+  std::string message;
+  return sparse_cholesky_->Factorize(lhs.get(), &message);
 }
 
 void VisibilityBasedPreconditioner::RightMultiply(const double* x,
                                                   double* y) const {
   CHECK_NOTNULL(x);
   CHECK_NOTNULL(y);
-  SuiteSparse* ss = const_cast<SuiteSparse*>(&ss_);
-
-  const int num_rows = m_->num_rows();
-  memcpy(CHECK_NOTNULL(tmp_rhs_)->x, x, m_->num_rows() * sizeof(*x));
-  // TODO(sameeragarwal): Better error handling.
-  std::string status;
-  cholmod_dense* solution =
-      CHECK_NOTNULL(ss->Solve(factor_, tmp_rhs_, &status));
-  memcpy(y, solution->x, sizeof(*y) * num_rows);
-  ss->Free(solution);
+  CHECK_NOTNULL(sparse_cholesky_.get());
+  std::string message;
+  sparse_cholesky_->Solve(x, y, &message);
 }
 
-int VisibilityBasedPreconditioner::num_rows() const {
-  return m_->num_rows();
-}
+int VisibilityBasedPreconditioner::num_rows() const { return m_->num_rows(); }
 
 // Classify camera/f_block pairs as in and out of the preconditioner,
 // based on whether the cluster pair that they belong to is in the
 // preconditioner or not.
 bool VisibilityBasedPreconditioner::IsBlockPairInPreconditioner(
-    const int block1,
-    const int block2) const {
+    const int block1, const int block2) const {
   int cluster1 = cluster_membership_[block1];
   int cluster2 = cluster_membership_[block2];
   if (cluster1 > cluster2) {
@@ -498,8 +454,7 @@ bool VisibilityBasedPreconditioner::IsBlockPairInPreconditioner(
 }
 
 bool VisibilityBasedPreconditioner::IsBlockPairOffDiagonal(
-    const int block1,
-    const int block2) const {
+    const int block1, const int block2) const {
   return (cluster_membership_[block1] != cluster_membership_[block2]);
 }
 
@@ -559,11 +514,13 @@ WeightedGraph<int>* VisibilityBasedPreconditioner::CreateClusterGraph(
 
   for (int i = 0; i < num_clusters_; ++i) {
     const set<int>& cluster_i = cluster_visibility[i];
-    for (int j = i+1; j < num_clusters_; ++j) {
+    for (int j = i + 1; j < num_clusters_; ++j) {
       vector<int> intersection;
       const set<int>& cluster_j = cluster_visibility[j];
-      set_intersection(cluster_i.begin(), cluster_i.end(),
-                       cluster_j.begin(), cluster_j.end(),
+      set_intersection(cluster_i.begin(),
+                       cluster_i.end(),
+                       cluster_j.begin(),
+                       cluster_j.end(),
                        back_inserter(intersection));
 
       if (intersection.size() > 0) {
@@ -614,9 +571,8 @@ void VisibilityBasedPreconditioner::FlattenMembershipMap(
       cluster_id = camera_id % num_clusters_;
     }
 
-    const int index = FindWithDefault(cluster_id_to_index,
-                                      cluster_id,
-                                      cluster_id_to_index.size());
+    const int index = FindWithDefault(
+        cluster_id_to_index, cluster_id, cluster_id_to_index.size());
 
     if (index == cluster_id_to_index.size()) {
       cluster_id_to_index[cluster_id] = index;
@@ -629,5 +585,3 @@ void VisibilityBasedPreconditioner::FlattenMembershipMap(
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/visibility_based_preconditioner.h

@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2017 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -49,15 +49,15 @@
 #define CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_
 
 #include <set>
-#include <vector>
 #include <utility>
+#include <vector>
 #include "ceres/collections_port.h"
 #include "ceres/graph.h"
 #include "ceres/internal/macros.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/linear_solver.h"
 #include "ceres/preconditioner.h"
-#include "ceres/suitesparse.h"
+#include "ceres/sparse_cholesky.h"
 
 namespace ceres {
 namespace internal {
@@ -122,8 +122,6 @@ class SchurEliminatorBase;
 //      *A.block_structure(), options);
 //   preconditioner.Update(A, NULL);
 //   preconditioner.RightMultiply(x, y);
-//
-#ifndef CERES_NO_SUITESPARSE
 class VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner {
  public:
   // Initialize the symbolic structure of the preconditioner. bs is
@@ -189,43 +187,9 @@ class VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner {
 
   // Preconditioner matrix.
   scoped_ptr<BlockRandomAccessSparseMatrix> m_;
-
-  // RightMultiply is a const method for LinearOperators. It is
-  // implemented using CHOLMOD's sparse triangular matrix solve
-  // function. This however requires non-const access to the
-  // SuiteSparse context object, even though it does not result in any
-  // of the state of the preconditioner being modified.
-  SuiteSparse ss_;
-
-  // Symbolic and numeric factorization of the preconditioner.
-  cholmod_factor* factor_;
-
-  // Temporary vector used by RightMultiply.
-  cholmod_dense* tmp_rhs_;
+  scoped_ptr<SparseCholesky> sparse_cholesky_;
   CERES_DISALLOW_COPY_AND_ASSIGN(VisibilityBasedPreconditioner);
 };
-#else  // SuiteSparse
-// If SuiteSparse is not compiled in, the preconditioner is not
-// available.
-class VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner {
- public:
-  VisibilityBasedPreconditioner(const CompressedRowBlockStructure& bs,
-                                const Preconditioner::Options& options) {
-    LOG(FATAL) << "Visibility based preconditioning is not available. Please "
-        "build Ceres with SuiteSparse.";
-  }
-  virtual ~VisibilityBasedPreconditioner() {}
-  virtual void RightMultiply(const double* x, double* y) const {}
-  virtual void LeftMultiply(const double* x, double* y) const {}
-  virtual int num_rows() const { return -1; }
-  virtual int num_cols() const { return -1; }
-
- private:
-  bool UpdateImpl(const BlockSparseMatrix& A, const double* D) {
-    return false;
-  }
-};
-#endif  // CERES_NO_SUITESPARSE
 
 }  // namespace internal
 }  // namespace ceres

internal/ceres/visibility_based_preconditioner_test.cc

@@ -28,11 +28,6 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/visibility_based_preconditioner.h"
 
 #include "Eigen/Dense"
@@ -345,5 +340,3 @@ namespace internal {
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

internal/ceres/visibility_test.cc

@@ -29,11 +29,6 @@
 // Author: kushalav@google.com (Avanish Kushal)
 //         sameeragarwal@google.com (Sameer Agarwal)
 
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
-#ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/visibility.h"
 
 #include <set>
@@ -209,5 +204,3 @@ TEST(VisibilityTest, NoEBlocks) {
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE

jni/Android.mk

@@ -145,6 +145,7 @@ LOCAL_SRC_FILES := $(CERES_SRC_PATH)/array_utils.cc \
                    $(CERES_SRC_PATH)/dynamic_compressed_row_jacobian_writer.cc \
                    $(CERES_SRC_PATH)/dynamic_compressed_row_sparse_matrix.cc \
                    $(CERES_SRC_PATH)/dynamic_sparse_normal_cholesky_solver.cc \
+                   $(CERES_SRC_PATH)/eigensparse.cc \
                    $(CERES_SRC_PATH)/evaluator.cc \
                    $(CERES_SRC_PATH)/file.cc \
                    $(CERES_SRC_PATH)/gradient_checker.cc \
@@ -186,6 +187,7 @@ LOCAL_SRC_FILES := $(CERES_SRC_PATH)/array_utils.cc \
                    $(CERES_SRC_PATH)/scratch_evaluate_preparer.cc \
                    $(CERES_SRC_PATH)/solver.cc \
                    $(CERES_SRC_PATH)/solver_utils.cc \
+                   $(CERES_SRC_PATH)/sparse_cholesky.cc \
                    $(CERES_SRC_PATH)/sparse_matrix.cc \
                    $(CERES_SRC_PATH)/sparse_normal_cholesky_solver.cc \
                    $(CERES_SRC_PATH)/split.cc \