diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index b113de5..572e92a 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -63,6 +63,7 @@
     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 @@
     solver.cc
     solver_utils.cc
     sparse_matrix.cc
+    sparse_cholesky.cc
     sparse_normal_cholesky_solver.cc
     split.cc
     stringprintf.cc
@@ -339,6 +341,7 @@
   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.
diff --git a/internal/ceres/canonical_views_clustering.cc b/internal/ceres/canonical_views_clustering.cc
index b655b1e..b3e9c22 100644
--- a/internal/ceres/canonical_views_clustering.cc
+++ b/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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/canonical_views_clustering.h b/internal/ceres/canonical_views_clustering.h
index 6b0b38a..0847f48 100644
--- a/internal/ceres/canonical_views_clustering.h
+++ b/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 @@
 }  // namespace internal
 }  // namespace ceres
 
-#endif  // CERES_NO_SUITESPARSE
 #endif  // CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_
diff --git a/internal/ceres/canonical_views_clustering_test.cc b/internal/ceres/canonical_views_clustering_test.cc
index 7282a3f..0c15fc7 100644
--- a/internal/ceres/canonical_views_clustering_test.cc
+++ b/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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/cxsparse.cc b/internal/ceres/cxsparse.cc
index 60b5808..9e9c599 100644
--- a/internal/ceres/cxsparse.cc
+++ b/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 @@
 
 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 @@
   }
 }
 
+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;
+        reinterpret_cast<CS_ENTRY*>(cs_di_malloc(num_cols, sizeof(CS_ENTRY)));
+    scratch_size_ = num_cols;
   }
 
-  // Solve using Cholesky factorization
-  csn* numeric_factorization = cs_di_chol(A, symbolic_factorization);
-  if (numeric_factorization == NULL) {
-    LOG(WARNING) << "Cholesky factorization failed.";
-    return false;
-  }
-
-  // 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 @@
 
   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 @@
   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_factor->pinv = cs_pinv(&scalar_ordering[0], A->n);
+  cs* permuted_A = cs_symperm(A, symbolic_factor->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->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 @@
   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);
 }
 
-void CXSparse::Free(cs_dis* symbolic_factorization) {
-  cs_di_sfree(symbolic_factorization);
+CompressedRowSparseMatrix::StorageType CXSparseCholesky::StorageType() const {
+  return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
+}
+
+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
diff --git a/internal/ceres/cxsparse.h b/internal/ceres/cxsparse.h
index 26dd192..14b1e01 100644
--- a/internal/ceres/cxsparse.h
+++ b/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 @@
   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 @@
 
   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 @@
   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
 
diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
index f6acd2c..4411223 100644
--- a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
@@ -211,20 +211,10 @@
   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
diff --git a/internal/ceres/eigensparse.cc b/internal/ceres/eigensparse.cc
new file mode 100644
index 0000000..aa6b6a9
--- /dev/null
+++ b/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
diff --git a/internal/ceres/eigensparse.h b/internal/ceres/eigensparse.h
new file mode 100644
index 0000000..402bb5c
--- /dev/null
+++ b/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_
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index fb9332c..abd630f 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -69,6 +69,10 @@
   LINEAR_SOLVER_FATAL_ERROR
 };
 
+enum OrderingType {
+  NATURAL,
+  AMD
+};
 
 class LinearOperator;
 
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 26c0e89..fa921b4 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/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 @@
     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 @@
 
 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 @@
   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;
   }
 
+  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);
+  }
+
   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);
-    }
-  } 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);
+  summary.termination_type = sparse_cholesky_->FactorAndSolve(
+      lhs.get(), rhs(), solution, &summary.message);
   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.";
-  }
-
-  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.
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h
index 714dafc..74665f0 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/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 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 @@
 // 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 @@
   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);
 };
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc
index cf9b09b..8228efb 100644
--- a/internal/ceres/schur_complement_solver_test.cc
+++ b/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 @@
     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 @@
   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);
 }
diff --git a/internal/ceres/single_linkage_clustering.cc b/internal/ceres/single_linkage_clustering.cc
index 490fcd8..9e9342a 100644
--- a/internal/ceres/single_linkage_clustering.cc
+++ b/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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/single_linkage_clustering.h b/internal/ceres/single_linkage_clustering.h
index fb02f01..8d1f02b 100644
--- a/internal/ceres/single_linkage_clustering.h
+++ b/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 @@
 }  // namespace internal
 }  // namespace ceres
 
-#endif  // CERES_NO_SUITESPARSE
 #endif  // CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_
diff --git a/internal/ceres/single_linkage_clustering_test.cc b/internal/ceres/single_linkage_clustering_test.cc
index 6bd4bd6..ca1a661 100644
--- a/internal/ceres/single_linkage_clustering_test.cc
+++ b/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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/sparse_cholesky.cc b/internal/ceres/sparse_cholesky.cc
new file mode 100644
index 0000000..1b5e638
--- /dev/null
+++ b/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
diff --git a/internal/ceres/sparse_cholesky.h b/internal/ceres/sparse_cholesky.h
new file mode 100644
index 0000000..7c2ba9d
--- /dev/null
+++ b/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_
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc
new file mode 100644
index 0000000..0c6e126
--- /dev/null
+++ b/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
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 557ebd8..3253548 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/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 @@
     }
     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
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index ef9d7be..597cff7 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/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 @@
       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);
 };
 
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index 1a0927b..024996c 100644
--- a/internal/ceres/suitesparse.cc
+++ b/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 @@
 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 @@
   return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_);
 }
 
-
 cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose(
     TripletSparseMatrix* A) {
   cholmod_triplet triplet;
@@ -127,12 +123,12 @@
 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 @@
   }
 
   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::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::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 @@
                                             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 @@
   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 @@
       *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::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 @@
 #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
 
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 380d76e..da02cc3 100644
--- a/internal/ceres/suitesparse.h
+++ b/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 @@
   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 @@
 
 typedef void cholmod_factor;
 
+namespace ceres {
+namespace internal {
+
 class SuiteSparse {
  public:
   // Defining this static function even when SuiteSparse is not
@@ -301,6 +325,9 @@
   void Free(void* arg) {}
 };
 
+}  // namespace internal
+}  // namespace ceres
+
 #endif  // CERES_NO_SUITESPARSE
 
 #endif  // CERES_INTERNAL_SUITESPARSE_H_
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc
index a670f00..401db1e 100644
--- a/internal/ceres/unsymmetric_linear_solver_test.cc
+++ b/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 @@
   options.dynamic_sparsity = true;
   TestSolver(options);
 }
-
 #endif  // CERES_USE_EIGEN_SPARSE
 
 }  // namespace internal
diff --git a/internal/ceres/visibility.cc b/internal/ceres/visibility.cc
index 55b5322..cb962a6 100644
--- a/internal/ceres/visibility.cc
+++ b/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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/visibility.h b/internal/ceres/visibility.h
index 142eb28..605682f 100644
--- a/internal/ceres/visibility.h
+++ b/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 @@
 }  // namespace internal
 }  // namespace ceres
 
-#endif  // CERES_NO_SUITESPARSE
 #endif  // CERES_INTERNAL_VISIBILITY_H_
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index 7e1ac67..90550b2 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/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 @@
 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 @@
   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 @@
   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 @@
     }
 
     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 @@
 // 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());
+  // S and convert it into a CompressedRowSparseMatrix.
+  const TripletSparseMatrix* tsm =
+      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);
+  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::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 @@
 
   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 @@
       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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h
index a627c13..40ce2c7 100644
--- a/internal/ceres/visibility_based_preconditioner.h
+++ b/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 @@
 //      *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 @@
 
   // 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
diff --git a/internal/ceres/visibility_based_preconditioner_test.cc b/internal/ceres/visibility_based_preconditioner_test.cc
index bd552fc..d2f13bc 100644
--- a/internal/ceres/visibility_based_preconditioner_test.cc
+++ b/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 ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/visibility_test.cc b/internal/ceres/visibility_test.cc
index a730b0e..73517da 100644
--- a/internal/ceres/visibility_test.cc
+++ b/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 @@
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_NO_SUITESPARSE
diff --git a/jni/Android.mk b/jni/Android.mk
index 9967885..6cffaae 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -145,6 +145,7 @@
                    $(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 @@
                    $(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 \