diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index d03c2ff..87e18e1 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -95,7 +95,7 @@
               "Options are: suite_sparse and cx_sparse.");
 DEFINE_string(dense_linear_algebra_library, "eigen",
               "Options are: eigen, lapack, and cuda");
-DEFINE_string(ordering, "automatic", "Options are: automatic, user.");
+DEFINE_string(ordering, "amd", "Options are: amd, nesdis and user");
 
 DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "
             "rotations. If false, angle axis is used.");
@@ -128,7 +128,6 @@
 DEFINE_string(initial_ply, "", "Export the BAL file data as a PLY file.");
 DEFINE_string(final_ply, "", "Export the refined BAL file data as a PLY "
               "file.");
-
 // clang-format on
 
 namespace ceres::examples {
@@ -228,24 +227,30 @@
   // the right ParameterBlock ordering, or by manually specifying a
   // suitable ordering vector and defining
   // Options::num_eliminate_blocks.
-  if (CERES_GET_FLAG(FLAGS_ordering) == "automatic") {
-    return;
+  if (CERES_GET_FLAG(FLAGS_ordering) == "user") {
+    auto* ordering = new ceres::ParameterBlockOrdering;
+
+    // The points come before the cameras.
+    for (int i = 0; i < num_points; ++i) {
+      ordering->AddElementToGroup(points + point_block_size * i, 0);
+    }
+
+    for (int i = 0; i < num_cameras; ++i) {
+      // When using axis-angle, there is a single parameter block for
+      // the entire camera.
+      ordering->AddElementToGroup(cameras + camera_block_size * i, 1);
+    }
+
+    options->linear_solver_ordering.reset(ordering);
+  } else {
+    if (CERES_GET_FLAG(FLAGS_ordering) == "amd") {
+      options->linear_solver_ordering_type = AMD;
+    } else if (CERES_GET_FLAG(FLAGS_ordering) == "nesdis") {
+      options->linear_solver_ordering_type = NESDIS;
+    } else {
+      LOG(FATAL) << "Unknown ordering type: " << CERES_GET_FLAG(FLAGS_ordering);
+    }
   }
-
-  auto* ordering = new ceres::ParameterBlockOrdering;
-
-  // The points come before the cameras.
-  for (int i = 0; i < num_points; ++i) {
-    ordering->AddElementToGroup(points + point_block_size * i, 0);
-  }
-
-  for (int i = 0; i < num_cameras; ++i) {
-    // When using axis-angle, there is a single parameter block for
-    // the entire camera.
-    ordering->AddElementToGroup(cameras + camera_block_size * i, 1);
-  }
-
-  options->linear_solver_ordering.reset(ordering);
 }
 
 void SetMinimizerOptions(Solver::Options* options) {
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 026fc1c0..43146b1 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -395,12 +395,33 @@
 #endif
 
     // The order in which variables are eliminated in a linear solver
-    // can have a significant of impact on the efficiency and accuracy
-    // of the method. e.g., when doing sparse Cholesky factorization,
+    // can have a significant impact on the efficiency and accuracy of
+    // the method. e.g., when doing sparse Cholesky factorization,
     // there are matrices for which a good ordering will give a
     // Cholesky factor with O(n) storage, where as a bad ordering will
     // result in an completely dense factor.
     //
+    // Sparse direct solvers like SPARSE_NORMAL_CHOLESKY and
+    // SPARSE_SCHUR use a fill reducing ordering of the columns and
+    // rows of the matrix being factorized before computing the
+    // numeric factorization.
+    //
+    // This enum controls the type of algorithm used to compute
+    // this fill reducing ordering. There is no single algorithm
+    // that works on all matrices, so determining which algorithm
+    // works better is a matter of empirical experimentation.
+    //
+    // Implementation status:
+    //
+    // AMD works for SUITE_SPARSE, CX_SPARSE, EIGEN_SPARSE &
+    // ACCELERATE_SPARSE.
+    //
+    // NESDIS currently works for SUITE_SPARSE when using
+    // SPARSE_NORMAL_CHOLESKY, SPARSE_SCHUR, CGNR + SUBSET,
+    // ITERATIVE_SCHUR + CLUSTER_JACOBI, ITERATIVE_SCHUR +
+    // CLUSTER_TRIDIAGONAL linear solvers.
+    LinearSolverOrderingType linear_solver_ordering_type = AMD;
+
     // Ceres allows the user to provide varying amounts of hints to
     // the solver about the variable elimination ordering to use. This
     // can range from no hints, where the solver is free to decide the
@@ -418,15 +439,47 @@
     // associated with it, that determines its order in the set of
     // groups.
     //
-    // Given such an ordering, Ceres ensures that the parameter blocks in
-    // the lowest numbered group are eliminated first, and then the
-    // parameter blocks in the next lowest numbered group and so on. Within
-    // each group, Ceres is free to order the parameter blocks as it
-    // chooses.
+    // The exact interpretation of this information depends on the
+    // values of linear_solver_ordering_type and
+    // linear_solver_type/preconditioner_type and
+    // sparse_linear_algebra_type.
     //
-    // If nullptr, then all parameter blocks are assumed to be in the
-    // same group and the solver is free to decide the best
-    // ordering.
+    // sparse_linear_algebra_library_type != SUITE_SPARSE
+    // ==================================================
+    //
+    // The value of linear_solver_ordering_type does not matter and:
+    //
+    // a. linear_solver_type = SPARSE_NORMAL_CHOLESKY or
+    //    linear_solver_type = CGNR and preconditioner_type = SUBSET
+    //
+    // then the value of linear_solver_ordering is ignored.
+    //
+    // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR
+    //
+    // then if linear_solver_ordering is non-null, then the lowest
+    // number group is used as the first elimination group to compute
+    // the Schur complement. All other groups are ignored.
+    //
+    // sparse_linear_algebra_library_type == SUITE_SPARSE
+    // ==================================================
+    //
+    // linear_solver_ordering_type = AMD
+    // ---------------------------------
+    //
+    // linear_solver_type = SPARSE_NORMAL_CHOLESKY or
+    // linear_solver_type = CGNR and preconditioner_type = SUBSET
+    //
+    // if linear_solver_ordering = nullptr, then Ceres assumes that
+    // all parameter blocks are in the same elimination group and uses
+    // the Approximate Minimum Degree algorithm to compute a fill
+    // reducing ordering.
+    //
+    // If linear_solver_order is not null, then a Constrained
+    // Approximate Minimum Degree (CAMD) ordering used where the
+    // parameter blocks in the lowest numbered group are eliminated
+    // first, and then the parameter blocks in the next lowest
+    // numbered group and so on. Within each group, CAMD free to order
+    // the parameter blocks as it chooses.
     //
     // e.g. Consider the linear system
     //
@@ -442,18 +495,53 @@
     //   {0: y}, {1: x} - eliminate y first.
     //   {0: x, y}      - Solver gets to decide the elimination order.
     //
-    // Thus, to have Ceres determine the ordering automatically using
-    // heuristics, put all the variables in group 0 and to control the
-    // ordering for every variable, create groups 0..N-1, one per
-    // variable, in the desired order.
+    // Thus, to have Ceres determine the ordering automatically, put
+    // all the variables in group 0 and to control the ordering for
+    // every variable, create groups 0..N-1, one per variable, in the
+    // desired order.
+    //
+    // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR
+    //
+    // If linear_solver_ordering is null then an greedy maximum
+    // independent set algorithm is used to compute the first
+    // elimination group, and AMD is applied to the corresponding
+    // Schur complement matrix or a subset thereof used that is used
+    // as a preconditioner.
+    //
+    // If linear_solver_ordering is not null, then CAMD is used to
+    // compute a fill reducing ordering.
+    //
+    // linear_solver_ordering = NESDIS
+    // -------------------------------
+    //
+    // a. linear_solver_type = SPARSE_NORMAL_CHOLESKY or
+    //    linear_solver_type = CGNR and preconditioner_type = SUBSET
+    //
+    // then the value of linear_solver_ordering is ignored and a
+    // Nested Dissection algorithm is used to compute a fill reducing
+    // ordering.
+    //
+    // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR
+    //
+    // If linear_solver_ordering is null then an greedy maximum
+    // independent set algorithm is used to compute the first
+    // elimination group, and Nested Dissection is applied to the
+    // corresponding Schur complement matrix or a subset thereof used
+    // that is used as a preconditioner.
+    //
+    // If linear_solver_ordering is not null, the parameter blocks in
+    // the lowest numbered group are eliminated to compute the Schur
+    // complement. All other groups are ignored and Nested Dissection
+    // is applied to the corresponding Schur complement matrix or a
+    // subset thereof used that is used as a preconditioner.
     //
     // Bundle Adjustment
     // -----------------
     //
-    // A particular case of interest is bundle adjustment, where the user
-    // has two options. The default is to not specify an ordering at all,
-    // the solver will see that the user wants to use a Schur type solver
-    // and figure out the right elimination ordering.
+    // A particular case of interest is bundle adjustment, where the
+    // user has two options. The default is not specifying ordering at
+    // all; in this case the solver will see that the user wants to
+    // use a Schur type solver and figure out an elimination ordering.
     //
     // But if the user already knows what parameter blocks are points and
     // what are cameras, they can save preprocessing time by partitioning
@@ -500,6 +588,10 @@
     // Jacobian matrix and generally speaking, there is no performance
     // penalty for doing so.
 
+    // TODO(sameeragarwal): Remove this option. It is too obscure and
+    // there is no clear way of figuring out when this is a useful
+    // thing to do.
+    //
     // In some rare cases, it is worth using a more complicated
     // reordering algorithm which has slightly better runtime
     // performance at the expense of an extra copy of the Jacobian
@@ -946,6 +1038,8 @@
         SPARSE_NORMAL_CHOLESKY;
 #endif
 
+    LinearSolverOrderingType linear_solver_ordering_type;
+
     // Size of the elimination groups given by the user as hints to
     // the linear solver.
     std::vector<int> linear_solver_ordering_given;
diff --git a/include/ceres/types.h b/include/ceres/types.h
index e522423..2fc5abc 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -183,6 +183,30 @@
   NO_SPARSE
 };
 
+// The order in which variables are eliminated in a linear solver
+// can have a significant of impact on the efficiency and accuracy
+// of the method. e.g., when doing sparse Cholesky factorization,
+// there are matrices for which a good ordering will give a
+// Cholesky factor with O(n) storage, where as a bad ordering will
+// result in an completely dense factor.
+//
+// So sparse direct solvers like SPARSE_NORMAL_CHOLESKY and
+// SPARSE_SCHUR and preconditioners like SUBSET, CLUSTER_JACOBI &
+// CLUSTER_TRIDIAGONAL use a fill reducing ordering of the columns and
+// rows of the matrix being factorized before actually the numeric
+// factorization.
+//
+// This enum controls the class of algorithm used to compute this
+// fill reducing ordering. There is no single algorithm that works
+// on all matrices, so determining which algorithm works better is a
+// matter of empirical experimentation.
+enum LinearSolverOrderingType {
+  // Approximate Minimum Degree.
+  AMD,
+  // Nested Dissection.
+  NESDIS
+};
+
 enum DenseLinearAlgebraLibraryType {
   EIGEN,
   LAPACK,
@@ -467,6 +491,11 @@
 CERES_EXPORT bool StringToSparseLinearAlgebraLibraryType(
     std::string value, SparseLinearAlgebraLibraryType* type);
 
+CERES_EXPORT const char* LinearSolverOrderingTypeToString(
+    LinearSolverOrderingType type);
+CERES_EXPORT bool StringToLinearSolverOrderingType(
+    std::string value, LinearSolverOrderingType* type);
+
 CERES_EXPORT const char* DenseLinearAlgebraLibraryTypeToString(
     DenseLinearAlgebraLibraryType type);
 CERES_EXPORT bool StringToDenseLinearAlgebraLibraryType(
diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc
index c6765fe..1f50a39 100644
--- a/internal/ceres/reorder_program.cc
+++ b/internal/ceres/reorder_program.cc
@@ -85,7 +85,6 @@
   return min_parameter_block_position;
 }
 
-#if defined(CERES_USE_EIGEN_SPARSE)
 Eigen::SparseMatrix<int> CreateBlockJacobian(
     const TripletSparseMatrix& block_jacobian_transpose) {
   using SparseMatrix = Eigen::SparseMatrix<int>;
@@ -105,9 +104,9 @@
   block_jacobian.setFromTriplets(triplets.begin(), triplets.end());
   return block_jacobian;
 }
-#endif
 
 void OrderingForSparseNormalCholeskyUsingSuiteSparse(
+    const LinearSolverOrderingType linear_solver_ordering_type,
     const TripletSparseMatrix& tsm_block_jacobian_transpose,
     const vector<ParameterBlock*>& parameter_blocks,
     const ParameterBlockOrdering& parameter_block_ordering,
@@ -120,40 +119,55 @@
   cholmod_sparse* block_jacobian_transpose = ss.CreateSparseMatrix(
       const_cast<TripletSparseMatrix*>(&tsm_block_jacobian_transpose));
 
-  // If the user did not supply a useful ordering, then just use
-  // regular AMD.
-  if (parameter_block_ordering.NumGroups() <= 1) {
-    ss.Ordering(block_jacobian_transpose, OrderingType::AMD, &ordering[0]);
-  } else {
-    vector<int> constraints;
-    for (auto* parameter_block : parameter_blocks) {
-      constraints.push_back(parameter_block_ordering.GroupId(
-          parameter_block->mutable_user_state()));
+  if (linear_solver_ordering_type == ceres::AMD) {
+    if (parameter_block_ordering.NumGroups() <= 1) {
+      // The user did not supply a useful ordering so just go ahead
+      // and use AMD.
+      ss.Ordering(block_jacobian_transpose, OrderingType::AMD, &ordering[0]);
+    } else {
+      // The user supplied an ordering, so use CAMD.
+      vector<int> constraints;
+      constraints.reserve(parameter_blocks.size());
+      for (auto* parameter_block : parameter_blocks) {
+        constraints.push_back(parameter_block_ordering.GroupId(
+            parameter_block->mutable_user_state()));
+      }
+
+      // Renumber the entries of constraints to be contiguous integers
+      // as CAMD requires that the group ids be in the range [0,
+      // parameter_blocks.size() - 1].
+      MapValuesToContiguousRange(constraints.size(), &constraints[0]);
+      ss.ConstrainedApproximateMinimumDegreeOrdering(
+          block_jacobian_transpose, &constraints[0], ordering);
     }
-
-    // Renumber the entries of constraints to be contiguous integers
-    // as CAMD requires that the group ids be in the range [0,
-    // parameter_blocks.size() - 1].
-    MapValuesToContiguousRange(constraints.size(), &constraints[0]);
-    ss.ConstrainedApproximateMinimumDegreeOrdering(
-        block_jacobian_transpose, &constraints[0], ordering);
+  } else if (linear_solver_ordering_type == ceres::NESDIS) {
+    // If nested dissection is chosen as an ordering algorithm, then
+    // ignore any user provided linear_solver_ordering.
+    CHECK(SuiteSparse::IsNestedDissectionAvailable())
+        << "Congratulations, you found a Ceres bug! "
+        << "Please report this error to the developers.";
+    ss.Ordering(block_jacobian_transpose, OrderingType::NESDIS, &ordering[0]);
+  } else {
+    LOG(FATAL) << "Congratulations, you found a Ceres bug! "
+               << "Please report this error to the developers.";
   }
 
-  VLOG(2) << "Block ordering stats: "
-          << " flops: " << ss.mutable_cc()->fl
-          << " lnz  : " << ss.mutable_cc()->lnz
-          << " anz  : " << ss.mutable_cc()->anz;
-
   ss.Free(block_jacobian_transpose);
 #endif  // CERES_NO_SUITESPARSE
 }
 
 void OrderingForSparseNormalCholeskyUsingCXSparse(
-    const TripletSparseMatrix& tsm_block_jacobian_transpose, int* ordering) {
+    const LinearSolverOrderingType linear_solver_ordering_type,
+    const TripletSparseMatrix& tsm_block_jacobian_transpose,
+    int* ordering) {
 #ifdef CERES_NO_CXSPARSE
   LOG(FATAL) << "Congratulations, you found a Ceres bug! "
              << "Please report this error to the developers.";
 #else
+  CHECK_NE(linear_solver_ordering_type, NESDIS)
+      << "Congratulations, you found a Ceres bug! "
+      << "Please report this error to the developers.";
+
   // CXSparse works with J'J instead of J'. So compute the block
   // sparsity for J'J and compute an approximate minimum degree
   // ordering.
@@ -173,13 +187,18 @@
 }
 
 void OrderingForSparseNormalCholeskyUsingEigenSparse(
-    const TripletSparseMatrix& tsm_block_jacobian_transpose, int* ordering) {
+    const LinearSolverOrderingType linear_solver_ordering_type,
+    const TripletSparseMatrix& tsm_block_jacobian_transpose,
+    int* ordering) {
 #ifndef CERES_USE_EIGEN_SPARSE
   LOG(FATAL) << "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.";
 #else
+  CHECK_NE(linear_solver_ordering_type, NESDIS)
+      << "Congratulations, you found a Ceres bug! "
+      << "Please report this error to the developers.";
 
   // This conversion from a TripletSparseMatrix to a Eigen::Triplet
   // matrix is unfortunate, but unavoidable for now. It is not a
@@ -326,7 +345,7 @@
 
 // Pre-order the columns corresponding to the Schur complement if
 // possible.
-static void MaybeReorderSchurComplementColumnsUsingSuiteSparse(
+static void ReorderSchurComplementColumnsUsingSuiteSparse(
     const ParameterBlockOrdering& parameter_block_ordering, Program* program) {
 #ifndef CERES_NO_SUITESPARSE
   SuiteSparse ss;
@@ -365,7 +384,7 @@
 #endif
 }
 
-static void MaybeReorderSchurComplementColumnsUsingEigen(
+static void ReorderSchurComplementColumnsUsingEigen(
     const int size_of_first_elimination_group,
     const ProblemImpl::ParameterMap& parameter_map,
     Program* program) {
@@ -420,6 +439,7 @@
 bool ReorderProgramForSchurTypeLinearSolver(
     const LinearSolverType linear_solver_type,
     const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const LinearSolverOrderingType linear_solver_ordering_type,
     const ProblemImpl::ParameterMap& parameter_map,
     ParameterBlockOrdering* parameter_block_ordering,
     Program* program,
@@ -485,12 +505,16 @@
   const int size_of_first_elimination_group =
       parameter_block_ordering->group_to_elements().begin()->second.size();
 
-  if (linear_solver_type == SPARSE_SCHUR) {
+  // Pre-ordering of the columns of the Schur complement only works if
+  // we are using approximate mininmum degree based ordering and
+  // SUITE_SPARSE or EIGEN_SPARSE.
+  if (linear_solver_type == SPARSE_SCHUR &&
+      linear_solver_ordering_type == ceres::AMD) {
     if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
-      MaybeReorderSchurComplementColumnsUsingSuiteSparse(
-          *parameter_block_ordering, program);
+      ReorderSchurComplementColumnsUsingSuiteSparse(*parameter_block_ordering,
+                                                    program);
     } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
-      MaybeReorderSchurComplementColumnsUsingEigen(
+      ReorderSchurComplementColumnsUsingEigen(
           size_of_first_elimination_group, parameter_map, program);
     }
   }
@@ -503,6 +527,7 @@
 
 bool ReorderProgramForSparseCholesky(
     const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const LinearSolverOrderingType linear_solver_ordering_type,
     const ParameterBlockOrdering& parameter_block_ordering,
     int start_row_block,
     Program* program,
@@ -526,12 +551,14 @@
 
   if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
     OrderingForSparseNormalCholeskyUsingSuiteSparse(
+        linear_solver_ordering_type,
         *tsm_block_jacobian_transpose,
         parameter_blocks,
         parameter_block_ordering,
         &ordering[0]);
   } else if (sparse_linear_algebra_library_type == CX_SPARSE) {
-    OrderingForSparseNormalCholeskyUsingCXSparse(*tsm_block_jacobian_transpose,
+    OrderingForSparseNormalCholeskyUsingCXSparse(linear_solver_ordering_type,
+                                                 *tsm_block_jacobian_transpose,
                                                  &ordering[0]);
   } else if (sparse_linear_algebra_library_type == ACCELERATE_SPARSE) {
     // Accelerate does not provide a function to perform reordering without
@@ -544,7 +571,9 @@
 
   } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
     OrderingForSparseNormalCholeskyUsingEigenSparse(
-        *tsm_block_jacobian_transpose, &ordering[0]);
+        linear_solver_ordering_type,
+        *tsm_block_jacobian_transpose,
+        &ordering[0]);
   }
 
   // Apply ordering.
@@ -569,4 +598,54 @@
   return it - residual_blocks->begin();
 }
 
+bool AreJacobianColumnsOrdered(
+    const LinearSolverType linear_solver_type,
+    const PreconditionerType preconditioner_type,
+    const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const LinearSolverOrderingType linear_solver_ordering_type,
+    const bool use_postordering,
+    const bool dynamic_sparsity) {
+  if (use_postordering || dynamic_sparsity) {
+    return false;
+  }
+
+  if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
+    if (linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
+      return true;
+    }
+    if (linear_solver_type == SPARSE_SCHUR &&
+        linear_solver_ordering_type == ceres::AMD) {
+      return true;
+    }
+  }
+
+  // For all sparse linear algebra libraries other than SuiteSparse,
+  // nested dissection is not used for pre-ordering.
+  if (linear_solver_ordering_type == ceres::NESDIS) {
+    return false;
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::EIGEN_SPARSE) {
+    if (linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        linear_solver_type == SPARSE_SCHUR ||
+        (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
+      return true;
+    }
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::CX_SPARSE) {
+    if (linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
+      return true;
+    }
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::ACCELERATE_SPARSE) {
+    return false;
+  }
+
+  return false;
+}
+
 }  // namespace ceres::internal
diff --git a/internal/ceres/reorder_program.h b/internal/ceres/reorder_program.h
index 64f4095..8ae456e 100644
--- a/internal/ceres/reorder_program.h
+++ b/internal/ceres/reorder_program.h
@@ -35,6 +35,7 @@
 
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
+#include "ceres/linear_solver.h"
 #include "ceres/parameter_block_ordering.h"
 #include "ceres/problem_impl.h"
 #include "ceres/types.h"
@@ -75,6 +76,7 @@
 CERES_NO_EXPORT bool ReorderProgramForSchurTypeLinearSolver(
     LinearSolverType linear_solver_type,
     SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    LinearSolverOrderingType linear_solver_ordering_type,
     const ProblemImpl::ParameterMap& parameter_map,
     ParameterBlockOrdering* parameter_block_ordering,
     Program* program,
@@ -92,6 +94,7 @@
 // ordering will take it into account, otherwise it will be ignored.
 CERES_NO_EXPORT bool ReorderProgramForSparseCholesky(
     SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    LinearSolverOrderingType linear_solver_ordering_type,
     const ParameterBlockOrdering& parameter_block_ordering,
     int start_row_block,
     Program* program,
@@ -111,6 +114,16 @@
     const std::unordered_set<ResidualBlockId>& bottom_residual_blocks,
     Program* program);
 
+// The return value of this function indicates whether the columns of
+// the Jacobian can be reordered using a fill reducing ordering.
+CERES_NO_EXPORT bool AreJacobianColumnsOrdered(
+    LinearSolverType linear_solver_type,
+    PreconditionerType preconditioner_type,
+    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    LinearSolverOrderingType linear_solver_ordering_type,
+    bool use_postordering,
+    bool dynamic_sparsity);
+
 }  // namespace ceres::internal
 
 #include "ceres/internal/reenable_warnings.h"
diff --git a/internal/ceres/reorder_program_test.cc b/internal/ceres/reorder_program_test.cc
index 0fbab4e..872f381 100644
--- a/internal/ceres/reorder_program_test.cc
+++ b/internal/ceres/reorder_program_test.cc
@@ -185,6 +185,7 @@
 
     std::string error;
     EXPECT_TRUE(ReorderProgramForSparseCholesky(ceres::SUITE_SPARSE,
+                                                ceres::AMD,
                                                 linear_solver_ordering,
                                                 0, /* use all rows */
                                                 program,
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index 150c555..dbaba3e 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -50,6 +50,7 @@
 #include "ceres/schur_templates.h"
 #include "ceres/solver_utils.h"
 #include "ceres/stringprintf.h"
+#include "ceres/suitesparse.h"
 #include "ceres/types.h"
 #include "ceres/wall_time.h"
 
@@ -107,6 +108,11 @@
   return true;
 }
 
+bool IsNestedDissectionAvailable(SparseLinearAlgebraLibraryType type) {
+  return (type == SUITE_SPARSE) &&
+         internal::SuiteSparse::IsNestedDissectionAvailable();
+}
+
 bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) {
   OPTION_GT(initial_trust_region_radius, 0.0);
   OPTION_GT(min_trust_region_radius, 0.0);
@@ -182,8 +188,10 @@
             name,
             sparse_linear_algebra_library_name);
         return false;
-      } else if (!IsSparseLinearAlgebraLibraryTypeAvailable(
-                     options.sparse_linear_algebra_library_type)) {
+      }
+
+      if (!IsSparseLinearAlgebraLibraryTypeAvailable(
+              options.sparse_linear_algebra_library_type)) {
         *error = StringPrintf(
             "Can't use %s with "
             "Solver::Options::sparse_linear_algebra_library_type = %s, "
@@ -192,6 +200,22 @@
             sparse_linear_algebra_library_name);
         return false;
       }
+
+      if (options.linear_solver_ordering_type == ceres::NESDIS &&
+          !IsNestedDissectionAvailable(
+              options.sparse_linear_algebra_library_type)) {
+        if (options.sparse_linear_algebra_library_type == SUITE_SPARSE) {
+          *error = StringPrintf(
+              "Can't use NESDIS with SUITE_SPARSE because SuiteSparse was "
+              "compiled without support for Metis.");
+        } else {
+          *error = StringPrintf(
+              "Can't use NESDIS with "
+              "Solver::Options::sparse_linear_algebra_library_type = %s.",
+              sparse_linear_algebra_library_name);
+        }
+        return false;
+      }
     }
   }
 
@@ -339,7 +363,7 @@
                                  &(summary->inner_iteration_ordering_given));
 
   // clang-format off
-  summary->dense_linear_algebra_library_type  = options.dense_linear_algebra_library_type;  //  NOLINT
+  summary->dense_linear_algebra_library_type  = options.dense_linear_algebra_library_type;  
   summary->dogleg_type                        = options.dogleg_type;
   summary->inner_iteration_time_in_seconds    = 0.0;
   summary->num_line_search_steps              = 0;
@@ -348,18 +372,19 @@
   summary->line_search_polynomial_minimization_time_in_seconds = 0.0;
   summary->line_search_total_time_in_seconds  = 0.0;
   summary->inner_iterations_given             = options.use_inner_iterations;
-  summary->line_search_direction_type         = options.line_search_direction_type;         //  NOLINT
-  summary->line_search_interpolation_type     = options.line_search_interpolation_type;     //  NOLINT
+  summary->line_search_direction_type         = options.line_search_direction_type;       
+  summary->line_search_interpolation_type     = options.line_search_interpolation_type;   
   summary->line_search_type                   = options.line_search_type;
   summary->linear_solver_type_given           = options.linear_solver_type;
   summary->max_lbfgs_rank                     = options.max_lbfgs_rank;
   summary->minimizer_type                     = options.minimizer_type;
-  summary->nonlinear_conjugate_gradient_type  = options.nonlinear_conjugate_gradient_type;  //  NOLINT
+  summary->nonlinear_conjugate_gradient_type  = options.nonlinear_conjugate_gradient_type; 
   summary->num_threads_given                  = options.num_threads;
   summary->preconditioner_type_given          = options.preconditioner_type;
-  summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; //  NOLINT
-  summary->trust_region_strategy_type         = options.trust_region_strategy_type;         //  NOLINT
-  summary->visibility_clustering_type         = options.visibility_clustering_type;         //  NOLINT
+  summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type;
+  summary->linear_solver_ordering_type        = options.linear_solver_ordering_type;       
+  summary->trust_region_strategy_type         = options.trust_region_strategy_type;        
+  summary->visibility_clustering_type         = options.visibility_clustering_type;
   // clang-format on
 }
 
@@ -367,11 +392,14 @@
                         Solver::Summary* summary) {
   internal::OrderingToGroupSizes(pp.options.linear_solver_ordering.get(),
                                  &(summary->linear_solver_ordering_used));
+  // TODO(sameeragarwal): Update the preprocessor to collapse the
+  // second and higher groups into one group when nested dissection is
+  // used.
   internal::OrderingToGroupSizes(pp.options.inner_iteration_ordering.get(),
                                  &(summary->inner_iteration_ordering_used));
 
   // clang-format off
-  summary->inner_iterations_used          = pp.inner_iteration_minimizer.get() != nullptr;     // NOLINT
+  summary->inner_iterations_used          = pp.inner_iteration_minimizer != nullptr;   
   summary->linear_solver_type_used        = pp.linear_solver_options.type;
   summary->num_threads_used               = pp.options.num_threads;
   summary->preconditioner_type_used       = pp.options.preconditioner_type;
@@ -379,7 +407,7 @@
 
   internal::SetSummaryFinalCost(summary);
 
-  if (pp.reduced_program.get() != nullptr) {
+  if (pp.reduced_program != nullptr) {
     SummarizeReducedProgram(*pp.reduced_program, summary);
   }
 
@@ -389,7 +417,7 @@
   // case if the preprocessor failed, or if the reduced problem did
   // not contain any parameter blocks. Thus, only extract the
   // evaluator statistics if one exists.
-  if (pp.evaluator.get() != nullptr) {
+  if (pp.evaluator != nullptr) {
     const map<string, CallStatistics>& evaluator_statistics =
         pp.evaluator->Statistics();
     {
@@ -411,7 +439,7 @@
   // Again, like the evaluator, there may or may not be a linear
   // solver from which we can extract run time statistics. In
   // particular the line search solver does not use a linear solver.
-  if (pp.linear_solver.get() != nullptr) {
+  if (pp.linear_solver != nullptr) {
     const map<string, CallStatistics>& linear_solver_statistics =
         pp.linear_solver->Statistics();
     const CallStatistics& call_stats = FindWithDefault(
@@ -666,17 +694,6 @@
                         dense_linear_algebra_library_type));
     }
 
-    if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
-        linear_solver_type_used == SPARSE_SCHUR ||
-        (linear_solver_type_used == ITERATIVE_SCHUR &&
-         (preconditioner_type_used == CLUSTER_JACOBI ||
-          preconditioner_type_used == CLUSTER_TRIDIAGONAL))) {
-      StringAppendF(&report,
-                    "\nSparse linear algebra library %15s\n",
-                    SparseLinearAlgebraLibraryTypeToString(
-                        sparse_linear_algebra_library_type));
-    }
-
     StringAppendF(&report,
                   "Trust region strategy     %19s",
                   TrustRegionStrategyTypeToString(trust_region_strategy_type));
@@ -687,9 +704,23 @@
         StringAppendF(&report, " (SUBSPACE)");
       }
     }
-    StringAppendF(&report, "\n");
-    StringAppendF(&report, "\n");
 
+    if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
+        linear_solver_type_used == SPARSE_SCHUR ||
+        (linear_solver_type_used == CGNR &&
+         preconditioner_type_used == SUBSET) ||
+        (linear_solver_type_used == ITERATIVE_SCHUR &&
+         (preconditioner_type_used == CLUSTER_JACOBI ||
+          preconditioner_type_used == CLUSTER_TRIDIAGONAL))) {
+      StringAppendF(
+          &report,
+          "\nSparse linear algebra library %15s + %s\n",
+          SparseLinearAlgebraLibraryTypeToString(
+              sparse_linear_algebra_library_type),
+          LinearSolverOrderingTypeToString(linear_solver_ordering_type));
+    }
+
+    StringAppendF(&report, "\n");
     StringAppendF(&report, "%45s    %21s\n", "Given", "Used");
     StringAppendF(&report,
                   "Linear solver       %25s%25s\n",
diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc
index 9467b5a..ed62aa9 100644
--- a/internal/ceres/trust_region_preprocessor.cc
+++ b/internal/ceres/trust_region_preprocessor.cc
@@ -113,6 +113,7 @@
     return ReorderProgramForSchurTypeLinearSolver(
         options.linear_solver_type,
         options.sparse_linear_algebra_library_type,
+        options.linear_solver_ordering_type,
         pp->problem->parameter_map(),
         options.linear_solver_ordering.get(),
         pp->reduced_program.get(),
@@ -123,6 +124,7 @@
       !options.dynamic_sparsity) {
     return ReorderProgramForSparseCholesky(
         options.sparse_linear_algebra_library_type,
+        options.linear_solver_ordering_type,
         *options.linear_solver_ordering,
         0, /* use all the rows of the jacobian */
         pp->reduced_program.get(),
@@ -138,6 +140,7 @@
 
     return ReorderProgramForSparseCholesky(
         options.sparse_linear_algebra_library_type,
+        options.linear_solver_ordering_type,
         *options.linear_solver_ordering,
         pp->linear_solver_options.subset_preconditioner_start_row_block,
         pp->reduced_program.get(),
@@ -200,6 +203,7 @@
       options.visibility_clustering_type;
   pp->linear_solver_options.sparse_linear_algebra_library_type =
       options.sparse_linear_algebra_library_type;
+
   pp->linear_solver_options.dense_linear_algebra_library_type =
       options.dense_linear_algebra_library_type;
   pp->linear_solver_options.use_explicit_schur_complement =
@@ -210,10 +214,6 @@
   pp->linear_solver_options.max_num_refinement_iterations =
       options.max_num_refinement_iterations;
   pp->linear_solver_options.num_threads = options.num_threads;
-  pp->linear_solver_options.ordering_type = OrderingType::NATURAL;
-  if (options.use_postordering) {
-    pp->linear_solver_options.ordering_type = OrderingType::AMD;
-  }
   pp->linear_solver_options.context = pp->problem->context();
 
   if (IsSchurType(pp->linear_solver_options.type)) {
@@ -227,20 +227,24 @@
     if (pp->linear_solver_options.elimination_groups.size() == 1) {
       pp->linear_solver_options.elimination_groups.push_back(0);
     }
+  }
 
-    if (options.linear_solver_type == SPARSE_SCHUR) {
-      // When using SPARSE_SCHUR, we ignore the user's postordering
-      // preferences if CX_SPARSE is the sparse linear algebra
-      // backend.
-      //
-      // This ensures that the linear solver does not assume that a
-      // fill-reducing pre-ordering has been done.
-      //
-      // TODO(sameeragarwal): Implement the reordering of parameter
-      // blocks for CX_SPARSE.
-      if (options.sparse_linear_algebra_library_type == CX_SPARSE) {
-        pp->linear_solver_options.ordering_type = OrderingType::AMD;
-      }
+  if (AreJacobianColumnsOrdered(options.linear_solver_type,
+                                options.preconditioner_type,
+                                options.sparse_linear_algebra_library_type,
+                                options.linear_solver_ordering_type,
+                                options.use_postordering,
+                                options.dynamic_sparsity)) {
+    pp->linear_solver_options.ordering_type = OrderingType::NATURAL;
+  } else {
+    if (options.linear_solver_ordering_type == ceres::AMD) {
+      pp->linear_solver_options.ordering_type = OrderingType::AMD;
+    } else if (options.linear_solver_ordering_type == ceres::NESDIS) {
+      pp->linear_solver_options.ordering_type = OrderingType::NESDIS;
+    } else {
+      LOG(FATAL) << "Congratulations you have found a bug in Ceres Solver."
+                 << " Please report this to the maintainers. : "
+                 << options.linear_solver_ordering_type;
     }
   }
 
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc
index 4824267..eea0312 100644
--- a/internal/ceres/types.cc
+++ b/internal/ceres/types.cc
@@ -124,6 +124,23 @@
   return false;
 }
 
+const char* LinearSolverOrderingTypeToString(LinearSolverOrderingType type) {
+  switch (type) {
+    CASESTR(AMD);
+    CASESTR(NESDIS);
+    default:
+      return "UNKNOWN";
+  }
+}
+
+bool StringToLinearSolverOrderingType(string value,
+                                      LinearSolverOrderingType* type) {
+  UpperCase(&value);
+  STRENUM(AMD);
+  STRENUM(NESDIS);
+  return false;
+}
+
 const char* DenseLinearAlgebraLibraryTypeToString(
     DenseLinearAlgebraLibraryType type) {
   switch (type) {
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index bde2814..0a7b19f 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -106,7 +106,7 @@
   LinearSolver::Options sparse_cholesky_options;
   sparse_cholesky_options.sparse_linear_algebra_library_type =
       options_.sparse_linear_algebra_library_type;
-  sparse_cholesky_options.ordering_type = options.ordering_type;
+  sparse_cholesky_options.ordering_type = options_.ordering_type;
   sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options);
 
   const time_t init_time = time(nullptr);