Nested dissection for ACCELERATE_SPARSE & EIGEN_SPARSE

Change-Id: Iec8ea6b0a537559b48b59bcfc91b94b58cb2070e
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 87e18e1..832e7e4 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -92,10 +92,11 @@
               "single_linkage, canonical_views");
 
 DEFINE_string(sparse_linear_algebra_library, "suite_sparse",
-              "Options are: suite_sparse and cx_sparse.");
+              "Options are: suite_sparse, cx_sparse, accelerate_sparse and eigen_sparse.");
 DEFINE_string(dense_linear_algebra_library, "eigen",
               "Options are: eigen, lapack, and cuda");
-DEFINE_string(ordering, "amd", "Options are: amd, nesdis and user");
+DEFINE_string(ordering_type, "amd", "Options are: amd, nesdis");
+DEFINE_string(linear_solver_ordering, "automatic", "Options are: automatic and user");
 
 DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "
             "rotations. If false, angle axis is used.");
@@ -147,6 +148,9 @@
   CHECK(StringToDenseLinearAlgebraLibraryType(
       CERES_GET_FLAG(FLAGS_dense_linear_algebra_library),
       &options->dense_linear_algebra_library_type));
+  CHECK(
+      StringToLinearSolverOrderingType(CERES_GET_FLAG(FLAGS_ordering_type),
+                                       &options->linear_solver_ordering_type));
   options->use_explicit_schur_complement =
       CERES_GET_FLAG(FLAGS_explicit_schur_complement);
   options->use_mixed_precision_solves =
@@ -222,12 +226,10 @@
   // ITERATIVE_SCHUR solvers make use of this specialized
   // structure.
   //
-  // This can either be done by specifying Options::ordering_type =
-  // ceres::SCHUR, in which case Ceres will automatically determine
-  // the right ParameterBlock ordering, or by manually specifying a
-  // suitable ordering vector and defining
-  // Options::num_eliminate_blocks.
-  if (CERES_GET_FLAG(FLAGS_ordering) == "user") {
+  // This can either be done by specifying a
+  // Options::linear_solver_ordering or having Ceres figure it out
+  // automatically using a greedy maximum independent set algorithm.
+  if (CERES_GET_FLAG(FLAGS_linear_solver_ordering) == "user") {
     auto* ordering = new ceres::ParameterBlockOrdering;
 
     // The points come before the cameras.
@@ -242,14 +244,6 @@
     }
 
     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);
-    }
   }
 }
 
diff --git a/internal/ceres/accelerate_sparse.cc b/internal/ceres/accelerate_sparse.cc
index a0d7e85..68a307b 100644
--- a/internal/ceres/accelerate_sparse.cc
+++ b/internal/ceres/accelerate_sparse.cc
@@ -139,8 +139,23 @@
 
 template <typename Scalar>
 typename AccelerateSparse<Scalar>::SymbolicFactorization
-AccelerateSparse<Scalar>::AnalyzeCholesky(ASSparseMatrix* A) {
-  return SparseFactor(SparseFactorizationCholesky, A->structure);
+AccelerateSparse<Scalar>::AnalyzeCholesky(OrderingType ordering_type,
+                                          ASSparseMatrix* A) {
+  SparseSymbolicFactorOptions sfoption;
+  sfoption.control = SparseDefaultControl;
+  sfoption.orderMethod = SparseOrderDefault;
+  sfoption.order = nullptr;
+  sfoption.ignoreRowsAndColumns = nullptr;
+  sfoption.malloc = malloc;
+  sfoption.free = free;
+  sfoption.reportError = nullptr;
+
+  if (ordering_type == OrderingType::AMD) {
+    sfoption.orderMethod = SparseOrderAMD;
+  } else if (ordering_type == OrderingType::NESDIS) {
+    sfoption.orderMethod = SparseOrderMetis;
+  }
+  return SparseFactor(SparseFactorizationCholesky, A->structure, sfoption);
 }
 
 template <typename Scalar>
@@ -207,7 +222,8 @@
   if (!symbolic_factor_) {
     symbolic_factor_ = std::make_unique<
         typename SparseTypesTrait<Scalar>::SymbolicFactorization>(
-        as_.AnalyzeCholesky(&as_lhs));
+        as_.AnalyzeCholesky(ordering_type_, &as_lhs));
+
     if (symbolic_factor_->status != SparseStatusOK) {
       *message = StringPrintf(
           "Apple Accelerate Failure : Symbolic factorisation failed: %s",
diff --git a/internal/ceres/accelerate_sparse.h b/internal/ceres/accelerate_sparse.h
index 29d78e8..81f5b94 100644
--- a/internal/ceres/accelerate_sparse.h
+++ b/internal/ceres/accelerate_sparse.h
@@ -91,7 +91,8 @@
   //       objects internally).
   ASSparseMatrix CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A);
   // Computes a symbolic factorisation of A that can be used in Solve().
-  SymbolicFactorization AnalyzeCholesky(ASSparseMatrix* A);
+  SymbolicFactorization AnalyzeCholesky(OrderingType ordering_type,
+                                        ASSparseMatrix* A);
   // Compute the numeric Cholesky factorization of A, given its
   // symbolic factorization.
   NumericFactorization Cholesky(ASSparseMatrix* A,
diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
index 2621690..1852969 100644
--- a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
+++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
@@ -116,6 +116,7 @@
 TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparseAMD) {
   TestSolver(SUITE_SPARSE, OrderingType::AMD);
 }
+
 #ifndef CERES_NO_METIS
 TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparseNESDIS) {
   TestSolver(SUITE_SPARSE, OrderingType::NESDIS);
@@ -130,9 +131,15 @@
 #endif
 
 #ifdef CERES_USE_EIGEN_SPARSE
-TEST_F(DynamicSparseNormalCholeskySolverTest, Eigen) {
+TEST_F(DynamicSparseNormalCholeskySolverTest, EigenAMD) {
   TestSolver(EIGEN_SPARSE, OrderingType::AMD);
 }
+
+#ifndef CERES_NO_METIS
+TEST_F(DynamicSparseNormalCholeskySolverTest, EigenNESDIS) {
+  TestSolver(EIGEN_SPARSE, OrderingType::NESDIS);
+}
+#endif
 #endif  // CERES_USE_EIGEN_SPARSE
 
 }  // namespace internal
diff --git a/internal/ceres/eigensparse.cc b/internal/ceres/eigensparse.cc
index 4c68c7a..7ef1210 100644
--- a/internal/ceres/eigensparse.cc
+++ b/internal/ceres/eigensparse.cc
@@ -34,8 +34,10 @@
 
 #ifdef CERES_USE_EIGEN_SPARSE
 
+#include <iostream>  // This is needed because MetisSupport depends on iostream.
 #include <sstream>
 
+#include "Eigen/MetisSupport"
 #include "Eigen/SparseCholesky"
 #include "Eigen/SparseCore"
 #include "ceres/compressed_row_sparse_matrix.h"
@@ -144,6 +146,9 @@
   using WithAMDOrdering = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
                                                 Eigen::Upper,
                                                 Eigen::AMDOrdering<int>>;
+  using WithMetisOrdering = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
+                                                  Eigen::Upper,
+                                                  Eigen::MetisOrdering<int>>;
   using WithNaturalOrdering =
       Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
                             Eigen::Upper,
@@ -151,9 +156,20 @@
 
   if (ordering_type == OrderingType::AMD) {
     return std::make_unique<EigenSparseCholeskyTemplate<WithAMDOrdering>>();
-  } else {
-    return std::make_unique<EigenSparseCholeskyTemplate<WithNaturalOrdering>>();
+#ifndef CERES_NO_METIS
+  } else if (ordering_type == OrderingType::NESDIS) {
+    return std::make_unique<EigenSparseCholeskyTemplate<WithMetisOrdering>>();
   }
+#else
+  } else {
+    LOG(FATAL)
+        << "Congratulations you have found a bug in Ceres Solver. Please "
+           "report it to the Ceres Solver developers.";
+    return nullptr;
+  }
+#endif
+
+  return std::make_unique<EigenSparseCholeskyTemplate<WithNaturalOrdering>>();
 }
 
 EigenSparseCholesky::~EigenSparseCholesky() = default;
@@ -163,15 +179,29 @@
   using WithAMDOrdering = Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>,
                                                 Eigen::Upper,
                                                 Eigen::AMDOrdering<int>>;
+  using WithMetisOrdering = Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>,
+                                                  Eigen::Upper,
+                                                  Eigen::MetisOrdering<int>>;
   using WithNaturalOrdering =
       Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>,
                             Eigen::Upper,
                             Eigen::NaturalOrdering<int>>;
   if (ordering_type == OrderingType::AMD) {
     return std::make_unique<EigenSparseCholeskyTemplate<WithAMDOrdering>>();
-  } else {
-    return std::make_unique<EigenSparseCholeskyTemplate<WithNaturalOrdering>>();
+#ifndef CERES_NO_METIS
+  } else if (ordering_type == OrderingType::NESDIS) {
+    return std::make_unique<EigenSparseCholeskyTemplate<WithMetisOrdering>>();
   }
+#else
+  } else {
+    LOG(FATAL)
+        << "Congratulations you have found a bug in Ceres Solver. Please "
+           "report it to the Ceres Solver developers.";
+    return nullptr;
+  }
+#endif
+
+  return std::make_unique<EigenSparseCholeskyTemplate<WithNaturalOrdering>>();
 }
 
 FloatEigenSparseCholesky::~FloatEigenSparseCholesky() = default;
diff --git a/internal/ceres/eigensparse.h b/internal/ceres/eigensparse.h
index e1a6b96..1ab1c57 100644
--- a/internal/ceres/eigensparse.h
+++ b/internal/ceres/eigensparse.h
@@ -48,6 +48,17 @@
 
 namespace ceres::internal {
 
+class EigenSparse {
+ public:
+  static constexpr bool IsNestedDissectionAvailable() noexcept {
+#ifdef CERES_NO_METIS
+    return false;
+#else
+    return true;
+#endif
+  }
+};
+
 class CERES_NO_EXPORT EigenSparseCholesky : public SparseCholesky {
  public:
   // Factory
diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc
index 2757725..b19a1b4 100644
--- a/internal/ceres/reorder_program.cc
+++ b/internal/ceres/reorder_program.cc
@@ -31,6 +31,7 @@
 #include "ceres/reorder_program.h"
 
 #include <algorithm>
+#include <iostream>  // Need this because MetisSupport refers to std::cerr.
 #include <memory>
 #include <numeric>
 #include <vector>
@@ -51,6 +52,7 @@
 #include "ceres/types.h"
 
 #ifdef CERES_USE_EIGEN_SPARSE
+#include "Eigen/MetisSupport"
 #include "Eigen/OrderingMethods"
 #endif
 
@@ -196,16 +198,13 @@
                 "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
-  // significant performance penalty in the grand scheme of
-  // things. The right thing to do here would be to get a compressed
-  // row sparse matrix representation of the jacobian and go from
-  // there. But that is a project for another day.
+  // TODO(sameeragarwal): This conversion from a TripletSparseMatrix
+  // to a Eigen::Triplet matrix is unfortunate, but unavoidable for
+  // now. It is not a significant performance penalty in the grand
+  // scheme of things. The right thing to do here would be to get a
+  // compressed row sparse matrix representation of the jacobian and
+  // go from there. But that is a project for another day.
   using SparseMatrix = Eigen::SparseMatrix<int>;
 
   const SparseMatrix block_jacobian =
@@ -213,9 +212,19 @@
   const SparseMatrix block_hessian =
       block_jacobian.transpose() * block_jacobian;
 
-  Eigen::AMDOrdering<int> amd_ordering;
   Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> perm;
-  amd_ordering(block_hessian, perm);
+  if (linear_solver_ordering_type == ceres::AMD) {
+    Eigen::AMDOrdering<int> amd_ordering;
+    amd_ordering(block_hessian, perm);
+  } else {
+#ifndef CERES_NO_METIS
+    perm.setIdentity(block_hessian.rows());
+#else
+    Eigen::MetisOrdering<int> metis_ordering;
+    metis_ordering(block_hessian, perm);
+#endif
+  }
+
   for (int i = 0; i < block_hessian.rows(); ++i) {
     ordering[i] = perm.indices()[i];
   }
@@ -385,13 +394,13 @@
 }
 
 static void ReorderSchurComplementColumnsUsingEigen(
+    LinearSolverOrderingType ordering_type,
     const int size_of_first_elimination_group,
     const ProblemImpl::ParameterMap& parameter_map,
     Program* program) {
 #if defined(CERES_USE_EIGEN_SPARSE)
   std::unique_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
       program->CreateJacobianBlockSparsityTranspose());
-
   using SparseMatrix = Eigen::SparseMatrix<int>;
   const SparseMatrix block_jacobian =
       CreateBlockJacobian(*tsm_block_jacobian_transpose);
@@ -412,9 +421,18 @@
   const SparseMatrix block_schur_complement =
       F.transpose() * F - F.transpose() * E * E.transpose() * F;
 
-  Eigen::AMDOrdering<int> amd_ordering;
   Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> perm;
-  amd_ordering(block_schur_complement, perm);
+  if (ordering_type == ceres::AMD) {
+    Eigen::AMDOrdering<int> amd_ordering;
+    amd_ordering(block_schur_complement, perm);
+  } else {
+#ifndef CERES_NO_METIS
+    perm.setIdentity(block_schur_complement.rows());
+#else
+    Eigen::MetisOrdering<int> metis_ordering;
+    metis_ordering(block_schur_complement, perm);
+#endif
+  }
 
   const vector<ParameterBlock*>& parameter_blocks = program->parameter_blocks();
   vector<ParameterBlock*> ordering(num_cols);
@@ -505,17 +523,21 @@
   const int size_of_first_elimination_group =
       parameter_block_ordering->group_to_elements().begin()->second.size();
 
-  // 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) {
+  if (linear_solver_type == SPARSE_SCHUR) {
+    if (sparse_linear_algebra_library_type == SUITE_SPARSE &&
+        linear_solver_ordering_type == ceres::AMD) {
+      // Preordering support for schur complement only works with AMD
+      // for now, since we are using CAMD.
+      //
+      // TODO(sameeragarwal): It maybe worth adding pre-ordering support for
+      // nested dissection too.
       ReorderSchurComplementColumnsUsingSuiteSparse(*parameter_block_ordering,
                                                     program);
     } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
-      ReorderSchurComplementColumnsUsingEigen(
-          size_of_first_elimination_group, parameter_map, program);
+      ReorderSchurComplementColumnsUsingEigen(linear_solver_ordering_type,
+                                              size_of_first_elimination_group,
+                                              parameter_map,
+                                              program);
     }
   }
 
@@ -612,11 +634,6 @@
         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;
   }
 
@@ -626,16 +643,24 @@
         (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
       return true;
     }
+    return false;
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::ACCELERATE_SPARSE) {
+    // Apple's accelerate framework does not allow direct access to
+    // ordering algorithms, so jacobian columns are never pre-ordered.
+    return false;
   }
 
   if (sparse_linear_algebra_library_type == ceres::CX_SPARSE) {
+    if (linear_solver_ordering_type == ceres::NESDIS) {
+      return false;
+    }
+
     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;
   }
 
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc
index 362a1f6..fe7a6b3 100644
--- a/internal/ceres/schur_complement_solver_test.cc
+++ b/internal/ceres/schur_complement_solver_test.cc
@@ -252,21 +252,37 @@
 #endif  // CERES_NO_CXSPARSE
 
 #ifndef CERES_NO_ACCELERATE_SPARSE
-// TODO(sameeragarwal): Extend these tests for NATURAL & NESDIS, once the linear
-// solver supports it.
-TEST_F(SchurComplementSolverTest, SparseSchurWithAccelerateSparseSmallProblem) {
+TEST_F(SchurComplementSolverTest,
+       SparseSchurWithAccelerateSparseSmallProblemAMD) {
   ComputeAndCompareSolutions(
       2, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
   ComputeAndCompareSolutions(
       2, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
 }
 
-TEST_F(SchurComplementSolverTest, SparseSchurWithAccelerateSparseLargeProblem) {
+TEST_F(SchurComplementSolverTest,
+       SparseSchurWithAccelerateSparseSmallProblemNESDIS) {
+  ComputeAndCompareSolutions(
+      2, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::NESDIS);
+  ComputeAndCompareSolutions(
+      2, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::NESDIS);
+}
+
+TEST_F(SchurComplementSolverTest,
+       SparseSchurWithAccelerateSparseLargeProblemAMD) {
   ComputeAndCompareSolutions(
       3, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
   ComputeAndCompareSolutions(
       3, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
 }
+
+TEST_F(SchurComplementSolverTest,
+       SparseSchurWithAccelerateSparseLargeProblemNESDIS) {
+  ComputeAndCompareSolutions(
+      3, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::NESDIS);
+  ComputeAndCompareSolutions(
+      3, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::NESDIS);
+}
 #endif  // CERES_NO_ACCELERATE_SPARSE
 
 #ifdef CERES_USE_EIGEN_SPARSE
@@ -277,6 +293,16 @@
       2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD);
 }
 
+#ifndef CERES_NO_METIS
+TEST_F(SchurComplementSolverTest,
+       SparseSchurWithEigenSparseSmallProblemNESDIS) {
+  ComputeAndCompareSolutions(
+      2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NESDIS);
+  ComputeAndCompareSolutions(
+      2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NESDIS);
+}
+#endif
+
 TEST_F(SchurComplementSolverTest,
        SparseSchurWithEigenSparseSmallProblemNATURAL) {
   ComputeAndCompareSolutions(
@@ -292,6 +318,16 @@
       3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD);
 }
 
+#ifndef CERES_NO_METIS
+TEST_F(SchurComplementSolverTest,
+       SparseSchurWithEigenSparseLargeProblemNESDIS) {
+  ComputeAndCompareSolutions(
+      3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NESDIS);
+  ComputeAndCompareSolutions(
+      3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NESDIS);
+}
+#endif
+
 TEST_F(SchurComplementSolverTest,
        SparseSchurWithEigenSparseLargeProblemNATURAL) {
   ComputeAndCompareSolutions(
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index dbaba3e..1f60fae 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -40,6 +40,7 @@
 #include "ceres/context.h"
 #include "ceres/context_impl.h"
 #include "ceres/detect_structure.h"
+#include "ceres/eigensparse.h"
 #include "ceres/gradient_checking_cost_function.h"
 #include "ceres/internal/export.h"
 #include "ceres/parameter_block_ordering.h"
@@ -109,8 +110,11 @@
 }
 
 bool IsNestedDissectionAvailable(SparseLinearAlgebraLibraryType type) {
-  return (type == SUITE_SPARSE) &&
-         internal::SuiteSparse::IsNestedDissectionAvailable();
+  return (((type == SUITE_SPARSE) &&
+           internal::SuiteSparse::IsNestedDissectionAvailable()) ||
+          (type == ACCELERATE_SPARSE) ||
+          ((type == EIGEN_SPARSE) &&
+           internal::EigenSparse::IsNestedDissectionAvailable()));
 }
 
 bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) {
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc
index 21b26a4..bc728e7 100644
--- a/internal/ceres/sparse_cholesky_test.cc
+++ b/internal/ceres/sparse_cholesky_test.cc
@@ -195,9 +195,9 @@
     SuiteSparseCholesky,
     SparseCholeskyTest,
     ::testing::Combine(::testing::Values(SUITE_SPARSE),
-                       ::testing::Values(OrderingType::NATURAL,
-                                         OrderingType::AMD,
-                                         OrderingType::NESDIS),
+                       ::testing::Values(OrderingType::AMD,
+                                         OrderingType::NESDIS,
+                                         OrderingType::NATURAL),
                        ::testing::Values(true, false)),
     ParamInfoToString);
 #endif
@@ -219,6 +219,7 @@
     SparseCholeskyTest,
     ::testing::Combine(::testing::Values(ACCELERATE_SPARSE),
                        ::testing::Values(OrderingType::AMD,
+                                         OrderingType::NESDIS,
                                          OrderingType::NATURAL),
                        ::testing::Values(true, false)),
     ParamInfoToString);
@@ -228,6 +229,7 @@
     SparseCholeskyTest,
     ::testing::Combine(::testing::Values(ACCELERATE_SPARSE),
                        ::testing::Values(OrderingType::AMD,
+                                         OrderingType::NESDIS,
                                          OrderingType::NATURAL),
                        ::testing::Values(true, false)),
     ParamInfoToString);
@@ -239,6 +241,7 @@
     SparseCholeskyTest,
     ::testing::Combine(::testing::Values(EIGEN_SPARSE),
                        ::testing::Values(OrderingType::AMD,
+                                         OrderingType::NESDIS,
                                          OrderingType::NATURAL),
                        ::testing::Values(true, false)),
     ParamInfoToString);
@@ -248,6 +251,7 @@
     SparseCholeskyTest,
     ::testing::Combine(::testing::Values(EIGEN_SPARSE),
                        ::testing::Values(OrderingType::AMD,
+                                         OrderingType::NESDIS,
                                          OrderingType::NATURAL),
                        ::testing::Values(true, false)),
     ParamInfoToString);