diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index 0d6f6bd..d93dd8d 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -29,9 +29,12 @@
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
 // This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
+#include "ceres/internal/config.h"
 
 #ifndef CERES_NO_SUITESPARSE
+
+#include <memory>
+#include <string>
 #include <vector>
 
 #include "ceres/compressed_col_sparse_matrix_utils.h"
@@ -41,11 +44,24 @@
 #include "ceres/triplet_sparse_matrix.h"
 #include "cholmod.h"
 
-namespace ceres {
-namespace internal {
+namespace ceres::internal {
+namespace {
+int OrderingTypeToCHOLMODEnum(OrderingType ordering_type) {
+  if (ordering_type == OrderingType::AMD) {
+    return CHOLMOD_AMD;
+  }
+  if (ordering_type == OrderingType::NESDIS) {
+    return CHOLMOD_NESDIS;
+  }
 
-using std::string;
-using std::vector;
+  if (ordering_type == OrderingType::NATURAL) {
+    return CHOLMOD_NATURAL;
+  }
+  LOG(FATAL) << "Congratulations you have discovered a bug in Ceres Solver."
+             << "Please report it to the developers. " << ordering_type;
+  return -1;
+}
+}  // namespace
 
 SuiteSparse::SuiteSparse() { cholmod_start(&cc_); }
 
@@ -102,9 +118,11 @@
   m.x = reinterpret_cast<void*>(A->mutable_values());
   m.z = nullptr;
 
-  if (A->storage_type() == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+  if (A->storage_type() ==
+      CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR) {
     m.stype = 1;
-  } else if (A->storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+  } else if (A->storage_type() ==
+             CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR) {
     m.stype = -1;
   } else {
     m.stype = 0;
@@ -143,19 +161,18 @@
 }
 
 cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A,
-                                             string* message) {
-  // Cholmod can try multiple re-ordering strategies to find a fill
-  // reducing ordering. Here we just tell it use AMD with automatic
-  // matrix dependence choice of supernodal versus simplicial
-  // factorization.
+                                             OrderingType ordering_type,
+                                             std::string* message) {
   cc_.nmethods = 1;
-  cc_.method[0].ordering = CHOLMOD_AMD;
-  cc_.supernodal = CHOLMOD_AUTO;
+  cc_.method[0].ordering = OrderingTypeToCHOLMODEnum(ordering_type);
+
+  // postordering with a NATURAL ordering leads to a significant regression in
+  // performance. See https://github.com/ceres-solver/ceres-solver/issues/905
+  if (ordering_type == OrderingType::NATURAL) {
+    cc_.postorder = 0;
+  }
 
   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 =
@@ -164,32 +181,22 @@
   }
 
   CHECK(factor != nullptr);
+  if (VLOG_IS_ON(2)) {
+    cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
+  }
+
   return factor;
 }
 
-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 nullptr;
-  }
-  return AnalyzeCholeskyWithUserOrdering(A, ordering, message);
-}
-
-cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(
-    cholmod_sparse* A, const vector<int>& ordering, string* message) {
+cholmod_factor* SuiteSparse::AnalyzeCholeskyWithGivenOrdering(
+    cholmod_sparse* A, const std::vector<int>& ordering, std::string* message) {
   CHECK_EQ(ordering.size(), A->nrow);
 
   cc_.nmethods = 1;
   cc_.method[0].ordering = CHOLMOD_GIVEN;
-
   cholmod_factor* factor =
-      cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), nullptr, 0, &cc_);
-  if (VLOG_IS_ON(2)) {
-    cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
-  }
+      cholmod_analyze_p(A, const_cast<int*>(ordering.data()), nullptr, 0, &cc_);
+
   if (cc_.status != CHOLMOD_OK) {
     *message =
         StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
@@ -197,40 +204,33 @@
   }
 
   CHECK(factor != nullptr);
-  return factor;
-}
-
-cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering(
-    cholmod_sparse* A, string* message) {
-  cc_.nmethods = 1;
-  cc_.method[0].ordering = CHOLMOD_NATURAL;
-  cc_.postorder = 0;
-
-  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);
-    return nullptr;
-  }
 
-  CHECK(factor != nullptr);
   return factor;
 }
 
-bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A,
-                                   const vector<int>& row_blocks,
-                                   const vector<int>& col_blocks,
-                                   vector<int>* ordering) {
+bool SuiteSparse::BlockOrdering(const cholmod_sparse* A,
+                                OrderingType ordering_type,
+                                const std::vector<Block>& row_blocks,
+                                const std::vector<Block>& col_blocks,
+                                std::vector<int>* ordering) {
+  if (ordering_type == OrderingType::NATURAL) {
+    ordering->resize(A->nrow);
+    for (int i = 0; i < A->nrow; ++i) {
+      (*ordering)[i] = i;
+    }
+    return true;
+  }
+
   const int num_row_blocks = row_blocks.size();
   const int num_col_blocks = col_blocks.size();
 
   // Arrays storing the compressed column structure of the matrix
-  // incoding the block sparsity of A.
-  vector<int> block_cols;
-  vector<int> block_rows;
+  // encoding the block sparsity of A.
+  std::vector<int> block_cols;
+  std::vector<int> block_rows;
 
   CompressedColumnScalarMatrixToBlockMatrix(reinterpret_cast<const int*>(A->i),
                                             reinterpret_cast<const int*>(A->p),
@@ -242,8 +242,8 @@
   block_matrix.nrow = num_row_blocks;
   block_matrix.ncol = num_col_blocks;
   block_matrix.nzmax = block_rows.size();
-  block_matrix.p = reinterpret_cast<void*>(&block_cols[0]);
-  block_matrix.i = reinterpret_cast<void*>(&block_rows[0]);
+  block_matrix.p = reinterpret_cast<void*>(block_cols.data());
+  block_matrix.i = reinterpret_cast<void*>(block_rows.data());
   block_matrix.x = nullptr;
   block_matrix.stype = A->stype;
   block_matrix.itype = CHOLMOD_INT;
@@ -252,8 +252,8 @@
   block_matrix.sorted = 1;
   block_matrix.packed = 1;
 
-  vector<int> block_ordering(num_row_blocks);
-  if (!cholmod_amd(&block_matrix, nullptr, 0, &block_ordering[0], &cc_)) {
+  std::vector<int> block_ordering(num_row_blocks);
+  if (!Ordering(&block_matrix, ordering_type, block_ordering.data())) {
     return false;
   }
 
@@ -261,9 +261,22 @@
   return true;
 }
 
+cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(
+    cholmod_sparse* A,
+    OrderingType ordering_type,
+    const std::vector<Block>& row_blocks,
+    const std::vector<Block>& col_blocks,
+    std::string* message) {
+  std::vector<int> ordering;
+  if (!BlockOrdering(A, ordering_type, row_blocks, col_blocks, &ordering)) {
+    return nullptr;
+  }
+  return AnalyzeCholeskyWithGivenOrdering(A, ordering, message);
+}
+
 LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A,
                                                   cholmod_factor* L,
-                                                  string* message) {
+                                                  std::string* message) {
   CHECK(A != nullptr);
   CHECK(L != nullptr);
 
@@ -281,48 +294,48 @@
   switch (cc_.status) {
     case CHOLMOD_NOT_INSTALLED:
       *message = "CHOLMOD failure: Method not installed.";
-      return LINEAR_SOLVER_FATAL_ERROR;
+      return LinearSolverTerminationType::FATAL_ERROR;
     case CHOLMOD_OUT_OF_MEMORY:
       *message = "CHOLMOD failure: Out of memory.";
-      return LINEAR_SOLVER_FATAL_ERROR;
+      return LinearSolverTerminationType::FATAL_ERROR;
     case CHOLMOD_TOO_LARGE:
       *message = "CHOLMOD failure: Integer overflow occurred.";
-      return LINEAR_SOLVER_FATAL_ERROR;
+      return LinearSolverTerminationType::FATAL_ERROR;
     case CHOLMOD_INVALID:
       *message = "CHOLMOD failure: Invalid input.";
-      return LINEAR_SOLVER_FATAL_ERROR;
+      return LinearSolverTerminationType::FATAL_ERROR;
     case CHOLMOD_NOT_POSDEF:
       *message = "CHOLMOD warning: Matrix not positive definite.";
-      return LINEAR_SOLVER_FAILURE;
+      return LinearSolverTerminationType::FAILURE;
     case CHOLMOD_DSMALL:
       *message =
           "CHOLMOD warning: D for LDL' or diag(L) or "
           "LL' has tiny absolute value.";
-      return LINEAR_SOLVER_FAILURE;
+      return LinearSolverTerminationType::FAILURE;
     case CHOLMOD_OK:
       if (cholmod_status != 0) {
-        return LINEAR_SOLVER_SUCCESS;
+        return LinearSolverTerminationType::SUCCESS;
       }
 
       *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;
+      return LinearSolverTerminationType::FATAL_ERROR;
     default:
       *message = StringPrintf(
           "Unknown cholmod return code: %d. "
           "Please report this to ceres-solver@googlegroups.com.",
           cc_.status);
-      return LINEAR_SOLVER_FATAL_ERROR;
+      return LinearSolverTerminationType::FATAL_ERROR;
   }
 
-  return LINEAR_SOLVER_FATAL_ERROR;
+  return LinearSolverTerminationType::FATAL_ERROR;
 }
 
 cholmod_dense* SuiteSparse::Solve(cholmod_factor* L,
                                   cholmod_dense* b,
-                                  string* message) {
+                                  std::string* message) {
   if (cc_.status != CHOLMOD_OK) {
     *message = "cholmod_solve failed. CHOLMOD status is not CHOLMOD_OK";
     return nullptr;
@@ -331,22 +344,34 @@
   return cholmod_solve(CHOLMOD_A, L, b, &cc_);
 }
 
-bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
-                                                   int* ordering) {
-  return cholmod_amd(matrix, nullptr, 0, ordering, &cc_);
+bool SuiteSparse::Ordering(cholmod_sparse* matrix,
+                           OrderingType ordering_type,
+                           int* ordering) {
+  CHECK_NE(ordering_type, OrderingType::NATURAL);
+  if (ordering_type == OrderingType::AMD) {
+    return cholmod_amd(matrix, nullptr, 0, ordering, &cc_);
+  }
+
+#ifdef CERES_NO_CHOLMOD_PARTITION
+  return false;
+#else
+  std::vector<int> CParent(matrix->nrow, 0);
+  std::vector<int> CMember(matrix->nrow, 0);
+  return cholmod_nested_dissection(
+      matrix, nullptr, 0, ordering, CParent.data(), CMember.data(), &cc_);
+#endif
 }
 
 bool SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering(
     cholmod_sparse* matrix, int* constraints, int* ordering) {
-#ifndef CERES_NO_CAMD
   return cholmod_camd(matrix, nullptr, 0, constraints, ordering, &cc_);
-#else
-  LOG(FATAL) << "Congratulations you have found a bug in Ceres."
-             << "Ceres Solver was compiled with SuiteSparse "
-             << "version 4.1.0 or less. Calling this function "
-             << "in that case is a bug. Please contact the"
-             << "the Ceres Solver developers.";
+}
+
+bool SuiteSparse::IsNestedDissectionAvailable() {
+#ifdef CERES_NO_CHOLMOD_PARTITION
   return false;
+#else
+  return true;
 #endif
 }
 
@@ -366,48 +391,61 @@
 }
 
 LinearSolverTerminationType SuiteSparseCholesky::Factorize(
-    CompressedRowSparseMatrix* lhs, string* message) {
+    CompressedRowSparseMatrix* lhs, std::string* message) {
   if (lhs == nullptr) {
-    *message = "Failure: Input lhs is NULL.";
-    return LINEAR_SOLVER_FATAL_ERROR;
+    *message = "Failure: Input lhs is nullptr.";
+    return LinearSolverTerminationType::FATAL_ERROR;
   }
 
   cholmod_sparse cholmod_lhs = ss_.CreateSparseMatrixTransposeView(lhs);
 
+  // If a factorization does not exist, compute the symbolic
+  // factorization first.
+  //
+  // If the ordering type is NATURAL, then there is no fill reducing
+  // ordering to be computed, regardless of block structure, so we can
+  // just call the scalar version of symbolic factorization. For
+  // SuiteSparse this is the common case since we have already
+  // pre-ordered the columns of the Jacobian.
+  //
+  // Similarly regardless of ordering type, if there is no block
+  // structure in the matrix we call the scalar version of symbolic
+  // factorization.
   if (factor_ == nullptr) {
-    if (ordering_type_ == NATURAL) {
-      factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&cholmod_lhs, message);
+    if (ordering_type_ == OrderingType::NATURAL ||
+        (lhs->col_blocks().empty() || lhs->row_blocks().empty())) {
+      factor_ = ss_.AnalyzeCholesky(&cholmod_lhs, ordering_type_, 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_ == nullptr) {
-      return LINEAR_SOLVER_FATAL_ERROR;
+      factor_ = ss_.BlockAnalyzeCholesky(&cholmod_lhs,
+                                         ordering_type_,
+                                         lhs->col_blocks(),
+                                         lhs->row_blocks(),
+                                         message);
     }
   }
 
+  if (factor_ == nullptr) {
+    return LinearSolverTerminationType::FATAL_ERROR;
+  }
+
+  // Compute and return the numeric factorization.
   return ss_.Cholesky(&cholmod_lhs, factor_, message);
 }
 
 CompressedRowSparseMatrix::StorageType SuiteSparseCholesky::StorageType()
     const {
-  return ((ordering_type_ == NATURAL)
-              ? CompressedRowSparseMatrix::UPPER_TRIANGULAR
-              : CompressedRowSparseMatrix::LOWER_TRIANGULAR);
+  return ((ordering_type_ == OrderingType::NATURAL)
+              ? CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR
+              : CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
 }
 
 LinearSolverTerminationType SuiteSparseCholesky::Solve(const double* rhs,
                                                        double* solution,
-                                                       string* message) {
+                                                       std::string* message) {
   // Error checking
   if (factor_ == nullptr) {
     *message = "Solve called without a call to Factorize first.";
-    return LINEAR_SOLVER_FATAL_ERROR;
+    return LinearSolverTerminationType::FATAL_ERROR;
   }
 
   const int num_cols = factor_->n;
@@ -416,15 +454,14 @@
       ss_.Solve(factor_, &cholmod_rhs, message);
 
   if (cholmod_dense_solution == nullptr) {
-    return LINEAR_SOLVER_FAILURE;
+    return LinearSolverTerminationType::FAILURE;
   }
 
   memcpy(solution, cholmod_dense_solution->x, num_cols * sizeof(*solution));
   ss_.Free(cholmod_dense_solution);
-  return LINEAR_SOLVER_SUCCESS;
+  return LinearSolverTerminationType::SUCCESS;
 }
 
-}  // namespace internal
-}  // namespace ceres
+}  // namespace ceres::internal
 
 #endif  // CERES_NO_SUITESPARSE