Refactor suitesparse.h/cc 1. Generalize SuiteSparse::AnalyzeCholesky and SuiteSparse::BlockAnalyzeCholesky from just doing AMD to taking OrderingType as an argument and using that to determine whether AMD & Nested Dissection algorithms are used for computing the fill-reducing ordering or a natural ordering when computing the symbolic factorization. 2. Remove AnalyzeCholeskyWithNaturalOrdering. 3. Replace and generalize SuiteSparse::BlockAMDOrdering with SuiteSparse::BlockOrdering which also takes OrderingType as an argument. Same for SuiteSparse::ApproximateMinimumDegreeOrdering and SuiteSparse::NestedDissectionOrdering by SuiteSparse::Ordering. 4. Remove LinearSolver::Options::use_postordering and replace it with LinearSolver::Options::ordering_type. 5. Replace Preconditioner::Options::use_postordering and replace it with Preconditioner::Options::ordering_type. 6. Add NESDIS to OrderingType. With the above changes, the linear solvers can now use Nested Dissection once this information is piped through the nonlinear solver. Change-Id: Ib8e93fbf34ae2981bf2ac54dcda9e25c7c213790
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc index c7cf4f8..12e2ef9 100644 --- a/internal/ceres/cgnr_solver.cc +++ b/internal/ceres/cgnr_solver.cc
@@ -80,7 +80,7 @@ options_.subset_preconditioner_start_row_block; preconditioner_options.sparse_linear_algebra_library_type = options_.sparse_linear_algebra_library_type; - preconditioner_options.use_postordering = options_.use_postordering; + preconditioner_options.ordering_type = options_.ordering_type; preconditioner_options.num_threads = options_.num_threads; preconditioner_options.context = options_.context; preconditioner_ =
diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc index f31004b..3abfb29 100644 --- a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc +++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
@@ -254,7 +254,8 @@ const int num_cols = A->num_cols(); cholmod_sparse lhs = ss.CreateSparseMatrixTransposeView(A); event_logger.AddEvent("Setup"); - cholmod_factor* factor = ss.AnalyzeCholesky(&lhs, &summary.message); + cholmod_factor* factor = + ss.AnalyzeCholesky(&lhs, options_.ordering_type, &summary.message); event_logger.AddEvent("Analysis"); if (factor == nullptr) {
diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc index f3bf650..2621690 100644 --- a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc +++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
@@ -93,12 +93,14 @@ } void TestSolver( - const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type) { + const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, + const OrderingType ordering_type) { LinearSolver::Options options; options.type = SPARSE_NORMAL_CHOLESKY; options.dynamic_sparsity = true; options.sparse_linear_algebra_library_type = sparse_linear_algebra_library_type; + options.ordering_type = ordering_type; ContextImpl context; options.context = &context; TestSolver(options, nullptr); @@ -111,20 +113,25 @@ }; #ifndef CERES_NO_SUITESPARSE -TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparse) { - TestSolver(SUITE_SPARSE); +TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparseAMD) { + TestSolver(SUITE_SPARSE, OrderingType::AMD); } +#ifndef CERES_NO_METIS +TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparseNESDIS) { + TestSolver(SUITE_SPARSE, OrderingType::NESDIS); +} +#endif #endif #ifndef CERES_NO_CXSPARSE TEST_F(DynamicSparseNormalCholeskySolverTest, CXSparse) { - TestSolver(CX_SPARSE); + TestSolver(CX_SPARSE, OrderingType::AMD); } #endif #ifdef CERES_USE_EIGEN_SPARSE TEST_F(DynamicSparseNormalCholeskySolverTest, Eigen) { - TestSolver(EIGEN_SPARSE); + TestSolver(EIGEN_SPARSE, OrderingType::AMD); } #endif // CERES_USE_EIGEN_SPARSE
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index 57c1a79..4916e41 100644 --- a/internal/ceres/linear_solver.h +++ b/internal/ceres/linear_solver.h
@@ -72,6 +72,7 @@ // the wrong arguments. FATAL_ERROR }; + inline std::ostream& operator<<(std::ostream& s, LinearSolverTerminationType type) { switch (type) { @@ -105,8 +106,27 @@ AMD, // Use the Approximate Minimum Degree algorithm to re-order // the matrix. + + NESDIS, // Use the Nested Dissection algorithm to re-order the matrix. }; +inline std::ostream& operator<<(std::ostream& s, OrderingType type) { + switch (type) { + case OrderingType::NATURAL: + s << "NATURAL"; + break; + case OrderingType::AMD: + s << "AMD"; + break; + case OrderingType::NESDIS: + s << "NESDIS"; + break; + default: + s << "UNKNOWN OrderingType"; + } + return s; +} + class LinearOperator; // Abstract base class for objects that implement algorithms for @@ -134,9 +154,9 @@ DenseLinearAlgebraLibraryType dense_linear_algebra_library_type = EIGEN; SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = SUITE_SPARSE; + OrderingType ordering_type = OrderingType::NATURAL; // See solver.h for information about these flags. - bool use_postordering = false; bool dynamic_sparsity = false; bool use_explicit_schur_complement = false;
diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h index fdc4d86..68b575f 100644 --- a/internal/ceres/preconditioner.h +++ b/internal/ceres/preconditioner.h
@@ -39,6 +39,7 @@ #include "ceres/internal/disable_warnings.h" #include "ceres/internal/export.h" #include "ceres/linear_operator.h" +#include "ceres/linear_solver.h" #include "ceres/sparse_matrix.h" #include "ceres/types.h" @@ -54,6 +55,7 @@ VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS; SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = SUITE_SPARSE; + OrderingType ordering_type = OrderingType::NATURAL; // When using the subset preconditioner, all row blocks starting // from this row block are used to construct the preconditioner. @@ -67,9 +69,6 @@ // and the preconditioner is the inverse of the matrix Q'Q. int subset_preconditioner_start_row_block = -1; - // See solver.h for information about these flags. - bool use_postordering = false; - // If possible, how many threads the preconditioner can use. int num_threads = 1;
diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc index bda4545..c6765fe 100644 --- a/internal/ceres/reorder_program.cc +++ b/internal/ceres/reorder_program.cc
@@ -123,7 +123,7 @@ // If the user did not supply a useful ordering, then just use // regular AMD. if (parameter_block_ordering.NumGroups() <= 1) { - ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose, &ordering[0]); + ss.Ordering(block_jacobian_transpose, OrderingType::AMD, &ordering[0]); } else { vector<int> constraints; for (auto* parameter_block : parameter_blocks) {
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc index 697d54a..362a1f6 100644 --- a/internal/ceres/schur_complement_solver_test.cc +++ b/internal/ceres/schur_complement_solver_test.cc
@@ -93,7 +93,7 @@ ceres::LinearSolverType linear_solver_type, ceres::DenseLinearAlgebraLibraryType dense_linear_algebra_library_type, ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, - bool use_postordering) { + ceres::internal::OrderingType ordering_type) { SetUpFromProblemId(problem_id); LinearSolver::Options options; options.elimination_groups.push_back(num_eliminate_blocks); @@ -104,7 +104,7 @@ dense_linear_algebra_library_type; options.sparse_linear_algebra_library_type = sparse_linear_algebra_library_type; - options.use_postordering = use_postordering; + options.ordering_type = ordering_type; ContextImpl context; options.context = &context; DetectStructure(*A->block_structure(), @@ -150,96 +150,154 @@ // 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); + ComputeAndCompareSolutions( + 2, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 2, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); } TEST_F(SchurComplementSolverTest, DenseSchurWithEigenLargeProblem) { - ComputeAndCompareSolutions(3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true); - ComputeAndCompareSolutions(3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true); + ComputeAndCompareSolutions( + 3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); } TEST_F(SchurComplementSolverTest, DenseSchurWithEigenVaryingFBlockSize) { - ComputeAndCompareSolutions(4, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true); + ComputeAndCompareSolutions( + 4, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); } #ifndef CERES_NO_LAPACK TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKSmallProblem) { - ComputeAndCompareSolutions(2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true); - ComputeAndCompareSolutions(2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true); + ComputeAndCompareSolutions( + 2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL); } TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKLargeProblem) { - ComputeAndCompareSolutions(3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true); - ComputeAndCompareSolutions(3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true); + ComputeAndCompareSolutions( + 3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL); } #endif #ifndef CERES_NO_SUITESPARSE TEST_F(SchurComplementSolverTest, - SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) { + SparseSchurWithSuiteSparseSmallProblemNATURAL) { ComputeAndCompareSolutions( - 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); - ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); -} - -TEST_F(SchurComplementSolverTest, - SparseSchurWithSuiteSparseSmallProblemPostOrdering) { - ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true); - ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true); -} - -TEST_F(SchurComplementSolverTest, - SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) { + 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); ComputeAndCompareSolutions( - 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); - ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); + 2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); } TEST_F(SchurComplementSolverTest, - SparseSchurWithSuiteSparseLargeProblemPostOrdering) { - ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true); - ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true); + SparseSchurWithSuiteSparseLargeProblemNATURAL) { + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL); } + +TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseSmallProblemAMD) { + ComputeAndCompareSolutions( + 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD); + ComputeAndCompareSolutions( + 2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD); +} + +TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseLargeProblemAMD) { + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD); + ComputeAndCompareSolutions( + 3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD); +} + +#ifndef CERES_NO_METIS +TEST_F(SchurComplementSolverTest, + SparseSchurWithSuiteSparseSmallProblemNESDIS) { + ComputeAndCompareSolutions( + 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS); + ComputeAndCompareSolutions( + 2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS); +} +TEST_F(SchurComplementSolverTest, + SparseSchurWithSuiteSparseLargeProblemNESDIS) { + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS); + ComputeAndCompareSolutions( + 3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS); +} +#endif // CERES_NO_METIS #endif // CERES_NO_SUITESPARSE #ifndef CERES_NO_CXSPARSE +// TODO(sameeragarwal): Extend these tests for NATURAL & NESDIS, once the linear +// solver supports it. TEST_F(SchurComplementSolverTest, SparseSchurWithCXSparseSmallProblem) { - ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true); - ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true); + ComputeAndCompareSolutions( + 2, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD); + ComputeAndCompareSolutions( + 2, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD); } TEST_F(SchurComplementSolverTest, SparseSchurWithCXSparseLargeProblem) { - ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true); - ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true); + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD); + ComputeAndCompareSolutions( + 3, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD); } #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) { ComputeAndCompareSolutions( - 2, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true); + 2, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD); ComputeAndCompareSolutions( - 2, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true); + 2, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD); } TEST_F(SchurComplementSolverTest, SparseSchurWithAccelerateSparseLargeProblem) { ComputeAndCompareSolutions( - 3, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true); + 3, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD); ComputeAndCompareSolutions( - 3, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true); + 3, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD); } #endif // CERES_NO_ACCELERATE_SPARSE #ifdef CERES_USE_EIGEN_SPARSE -TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseSmallProblem) { - ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true); - ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true); +TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseSmallProblemAMD) { + ComputeAndCompareSolutions( + 2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD); + ComputeAndCompareSolutions( + 2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD); } -TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseLargeProblem) { - ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true); - ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true); +TEST_F(SchurComplementSolverTest, + SparseSchurWithEigenSparseSmallProblemNATURAL) { + ComputeAndCompareSolutions( + 2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL); +} + +TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseLargeProblemAMD) { + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD); + ComputeAndCompareSolutions( + 3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD); +} + +TEST_F(SchurComplementSolverTest, + SparseSchurWithEigenSparseLargeProblemNATURAL) { + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL); + ComputeAndCompareSolutions( + 3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL); } #endif // CERES_USE_EIGEN_SPARSE
diff --git a/internal/ceres/sparse_cholesky.cc b/internal/ceres/sparse_cholesky.cc index ac7c9b2..c68b833 100644 --- a/internal/ceres/sparse_cholesky.cc +++ b/internal/ceres/sparse_cholesky.cc
@@ -44,17 +44,16 @@ std::unique_ptr<SparseCholesky> SparseCholesky::Create( const LinearSolver::Options& options) { - const OrderingType ordering_type = - options.use_postordering ? OrderingType::AMD : OrderingType::NATURAL; std::unique_ptr<SparseCholesky> sparse_cholesky; switch (options.sparse_linear_algebra_library_type) { case SUITE_SPARSE: #ifndef CERES_NO_SUITESPARSE if (options.use_mixed_precision_solves) { - sparse_cholesky = FloatSuiteSparseCholesky::Create(ordering_type); + sparse_cholesky = + FloatSuiteSparseCholesky::Create(options.ordering_type); } else { - sparse_cholesky = SuiteSparseCholesky::Create(ordering_type); + sparse_cholesky = SuiteSparseCholesky::Create(options.ordering_type); } break; #else @@ -64,9 +63,10 @@ case EIGEN_SPARSE: #ifdef CERES_USE_EIGEN_SPARSE if (options.use_mixed_precision_solves) { - sparse_cholesky = FloatEigenSparseCholesky::Create(ordering_type); + sparse_cholesky = + FloatEigenSparseCholesky::Create(options.ordering_type); } else { - sparse_cholesky = EigenSparseCholesky::Create(ordering_type); + sparse_cholesky = EigenSparseCholesky::Create(options.ordering_type); } break; #else @@ -77,9 +77,9 @@ case CX_SPARSE: #ifndef CERES_NO_CXSPARSE if (options.use_mixed_precision_solves) { - sparse_cholesky = FloatCXSparseCholesky::Create(ordering_type); + sparse_cholesky = FloatCXSparseCholesky::Create(options.ordering_type); } else { - sparse_cholesky = CXSparseCholesky::Create(ordering_type); + sparse_cholesky = CXSparseCholesky::Create(options.ordering_type); } break; #else @@ -89,10 +89,11 @@ case ACCELERATE_SPARSE: #ifndef CERES_NO_ACCELERATE_SPARSE if (options.use_mixed_precision_solves) { - sparse_cholesky = AppleAccelerateCholesky<float>::Create(ordering_type); + sparse_cholesky = + AppleAccelerateCholesky<float>::Create(options.ordering_type); } else { sparse_cholesky = - AppleAccelerateCholesky<double>::Create(ordering_type); + AppleAccelerateCholesky<double>::Create(options.ordering_type); } break; #else
diff --git a/internal/ceres/sparse_cholesky.h b/internal/ceres/sparse_cholesky.h index f4f8996..feea7aa 100644 --- a/internal/ceres/sparse_cholesky.h +++ b/internal/ceres/sparse_cholesky.h
@@ -63,7 +63,7 @@ // CompressedRowSparseMatrix lhs = ...; // std::string message; // CHECK_EQ(sparse_cholesky->Factorize(&lhs, &message), -// LinearSolverTerminationType::SUCCESS); +// LinearSolverTerminationType::SUCCESS); // Vector rhs = ...; // Vector solution = ...; // CHECK_EQ(sparse_cholesky->Solve(rhs.data(), solution.data(), &message),
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc index c9c8363..21b26a4 100644 --- a/internal/ceres/sparse_cholesky_test.cc +++ b/internal/ceres/sparse_cholesky_test.cc
@@ -114,8 +114,7 @@ LinearSolver::Options sparse_cholesky_options; sparse_cholesky_options.sparse_linear_algebra_library_type = sparse_linear_algebra_library_type; - sparse_cholesky_options.use_postordering = - (ordering_type == OrderingType::AMD); + sparse_cholesky_options.ordering_type = ordering_type; auto sparse_cholesky = SparseCholesky::Create(sparse_cholesky_options); const CompressedRowSparseMatrix::StorageType storage_type = sparse_cholesky->StorageType(); @@ -156,8 +155,7 @@ Param param = info.param; std::stringstream ss; ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_" - << (::testing::get<1>(param) == OrderingType::AMD ? "AMD" : "NATURAL") - << "_" + << ::testing::get<1>(param) << "_" << (::testing::get<2>(param) ? "UseBlockStructure" : "NoBlockStructure"); return ss.str(); } @@ -197,8 +195,9 @@ SuiteSparseCholesky, SparseCholeskyTest, ::testing::Combine(::testing::Values(SUITE_SPARSE), - ::testing::Values(OrderingType::AMD, - OrderingType::NATURAL), + ::testing::Values(OrderingType::NATURAL, + OrderingType::AMD, + OrderingType::NESDIS), ::testing::Values(true, false)), ParamInfoToString); #endif
diff --git a/internal/ceres/sparse_normal_cholesky_solver_test.cc b/internal/ceres/sparse_normal_cholesky_solver_test.cc index 6b3316c..39d9d82 100644 --- a/internal/ceres/sparse_normal_cholesky_solver_test.cc +++ b/internal/ceres/sparse_normal_cholesky_solver_test.cc
@@ -111,7 +111,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = SUITE_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = false; + options.ordering_type = OrderingType::NATURAL; ContextImpl context; options.context = &context; TestSolver(options); @@ -122,7 +122,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = SUITE_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = true; + options.ordering_type = OrderingType::AMD; ContextImpl context; options.context = &context; TestSolver(options); @@ -135,7 +135,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = CX_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = false; + options.ordering_type = OrderingType::NATURAL; ContextImpl context; options.context = &context; TestSolver(options); @@ -146,7 +146,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = CX_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = true; + options.ordering_type = OrderingType::AMD; ContextImpl context; options.context = &context; TestSolver(options); @@ -159,7 +159,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = ACCELERATE_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = false; + options.ordering_type = OrderingType::NATURAL; ContextImpl context; options.context = &context; TestSolver(options); @@ -170,7 +170,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = ACCELERATE_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = true; + options.ordering_type = OrderingType::AMD; ContextImpl context; options.context = &context; TestSolver(options); @@ -183,7 +183,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = EIGEN_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = false; + options.ordering_type = OrderingType::NATURAL; ContextImpl context; options.context = &context; TestSolver(options); @@ -194,7 +194,7 @@ LinearSolver::Options options; options.sparse_linear_algebra_library_type = EIGEN_SPARSE; options.type = SPARSE_NORMAL_CHOLESKY; - options.use_postordering = true; + options.ordering_type = OrderingType::AMD; ContextImpl context; options.context = &context; TestSolver(options);
diff --git a/internal/ceres/subset_preconditioner.cc b/internal/ceres/subset_preconditioner.cc index c804274..a228545 100644 --- a/internal/ceres/subset_preconditioner.cc +++ b/internal/ceres/subset_preconditioner.cc
@@ -51,7 +51,7 @@ LinearSolver::Options sparse_cholesky_options; sparse_cholesky_options.sparse_linear_algebra_library_type = options_.sparse_linear_algebra_library_type; - sparse_cholesky_options.use_postordering = options_.use_postordering; + sparse_cholesky_options.ordering_type = options_.ordering_type; sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options); }
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc index a7d28e6..82d20d4 100644 --- a/internal/ceres/suitesparse.cc +++ b/internal/ceres/suitesparse.cc
@@ -43,6 +43,23 @@ #include "cholmod.h" 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; + } + + 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 using std::string; using std::vector; @@ -145,42 +162,38 @@ } cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A, + OrderingType ordering_type, 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. cc_.nmethods = 1; - cc_.method[0].ordering = CHOLMOD_AMD; + cc_.method[0].ordering = OrderingTypeToCHOLMODEnum(ordering_type); cc_.supernodal = CHOLMOD_AUTO; - cholmod_factor* factor = cholmod_analyze(A, &cc_); - if (VLOG_IS_ON(2)) { - cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); + // TODO(sameeragarwal): The following modification exists entirely + // to preserve existing behaviour. It is not clear if this is + // actually the optimal thing to do. Once I understand the use of + // postordering better, I will come back and document this/remove it + // as the case maybe. + if (ordering_type == OrderingType::NATURAL) { + cc_.postorder = 0; } + cholmod_factor* factor = cholmod_analyze(A, &cc_); if (cc_.status != CHOLMOD_OK) { *message = StringPrintf("cholmod_analyze failed. error code: %d", cc_.status); return nullptr; } - CHECK(factor != nullptr); - 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; + if (VLOG_IS_ON(2)) { + cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); } - return AnalyzeCholeskyWithUserOrdering(A, ordering, message); } -cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering( +return factor; +} + +cholmod_factor* SuiteSparse::AnalyzeCholeskyWithGivenOrdering( cholmod_sparse* A, const vector<int>& ordering, string* message) { CHECK_EQ(ordering.size(), A->nrow); @@ -189,9 +202,6 @@ 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_); - } if (cc_.status != CHOLMOD_OK) { *message = StringPrintf("cholmod_analyze failed. error code: %d", cc_.status); @@ -199,33 +209,26 @@ } 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 vector<int>& row_blocks, + const vector<int>& col_blocks, + 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(); @@ -255,7 +258,7 @@ block_matrix.packed = 1; vector<int> block_ordering(num_row_blocks); - if (!cholmod_amd(&block_matrix, nullptr, 0, &block_ordering[0], &cc_)) { + if (!Ordering(&block_matrix, ordering_type, block_ordering.data())) { return false; } @@ -263,6 +266,18 @@ return true; } +cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(cholmod_sparse* A, + OrderingType ordering_type, + const vector<int>& row_blocks, + const vector<int>& col_blocks, + string* message) { + 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) { @@ -333,13 +348,14 @@ 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_); + } -bool SuiteSparse::NestedDissectionOrdering(cholmod_sparse* matrix, - int* ordering) { #ifdef CERES_NO_METIS return false; #else @@ -355,6 +371,14 @@ return cholmod_camd(matrix, nullptr, 0, constraints, ordering, &cc_); } +bool SuiteSparse::IsNestedDissectionAvailable() { +#ifdef CERES_NO_METIS + return false; +#else + return true; +#endif +} + std::unique_ptr<SparseCholesky> SuiteSparseCholesky::Create( const OrderingType ordering_type) { return std::unique_ptr<SparseCholesky>( @@ -379,23 +403,36 @@ 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_ == OrderingType::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 LinearSolverTerminationType::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); }
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h index b8dc666..e0c485d 100644 --- a/internal/ceres/suitesparse.h +++ b/internal/ceres/suitesparse.h
@@ -119,12 +119,11 @@ cholmod_sdmult(A, 0, alpha_, beta_, x, y, &cc_); } - // Find an ordering of A or AA' (if A is unsymmetric) that minimizes - // the fill-in in the Cholesky factorization of the corresponding - // matrix. This is done by using the AMD algorithm. - // - // Using this ordering, the symbolic Cholesky factorization of A (or - // AA') is computed and returned. + // Compute a symbolic factorization for A or AA' (if A is + // unsymmetric). If ordering_type is NATURAL, then no fill reducing + // ordering is computed, otherwise depending on the value of + // ordering_type AMD or Nested Dissection is used to compute a fill + // reducing ordering before the symbolic factorization is computed. // // A is not modified, only the pattern of non-zeros of A is used, // the actual numerical values in A are of no consequence. @@ -132,9 +131,13 @@ // message contains an explanation of the failures if any. // // Caller owns the result. - cholmod_factor* AnalyzeCholesky(cholmod_sparse* A, std::string* message); + cholmod_factor* AnalyzeCholesky(cholmod_sparse* A, + OrderingType ordering_type, + std::string* message); + // Block oriented version of AnalyzeCholesky. cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A, + OrderingType ordering_type, const std::vector<int>& row_blocks, const std::vector<int>& col_blocks, std::string* message); @@ -150,20 +153,11 @@ // message contains an explanation of the failures if any. // // Caller owns the result. - cholmod_factor* AnalyzeCholeskyWithUserOrdering( + cholmod_factor* AnalyzeCholeskyWithGivenOrdering( cholmod_sparse* A, const std::vector<int>& ordering, std::string* message); - // Perform a symbolic factorization of A without re-ordering A. No - // postordering of the elimination tree is performed. This ensures - // that the symbolic factor does not introduce an extra permutation - // on the matrix. See the documentation for CHOLMOD for more details. - // - // message contains an explanation of the failures if any. - cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A, - std::string* message); - // Use the symbolic factorization in L, to find the numerical // factorization for the matrix A or AA^T. Return true if // successful, false otherwise. L contains the numeric factorization @@ -183,50 +177,39 @@ cholmod_dense* b, std::string* message); + // Find a fill reducing ordering. ordering is expected to be large + // enough to hold the ordering. ordering_type must be AMD or NESDIS. + bool Ordering(cholmod_sparse* matrix, + OrderingType ordering_type, + int* ordering); + + // Find the block oriented fill reducing ordering of a matrix A, + // whose row and column blocks are given by row_blocks, and + // col_blocks respectively. The matrix may or may not be + // symmetric. The entries of col_blocks do not need to sum to the + // number of columns in A. If this is the case, only the first + // sum(col_blocks) are used to compute the ordering. + // // By virtue of the modeling layer in Ceres being block oriented, // all the matrices used by Ceres are also block oriented. When // doing sparse direct factorization of these matrices the - // fill-reducing ordering algorithms (in particular AMD) can either - // be run on the block or the scalar form of these matrices. The two - // SuiteSparse::AnalyzeCholesky methods allows the client to - // compute the symbolic factorization of a matrix by either using - // AMD on the matrix or a user provided ordering of the rows. - // - // But since the underlying matrices are block oriented, it is worth - // running AMD on just the block structure of these matrices and then - // lifting these block orderings to a full scalar ordering. This - // preserves the block structure of the permuted matrix, and exposes - // more of the super-nodal structure of the matrix to the numerical - // factorization routines. - // - // Find the block oriented AMD ordering of a matrix A, whose row and - // column blocks are given by row_blocks, and col_blocks - // respectively. The matrix may or may not be symmetric. The entries - // of col_blocks do not need to sum to the number of columns in - // A. If this is the case, only the first sum(col_blocks) are used - // to compute the ordering. - bool BlockAMDOrdering(const cholmod_sparse* A, - const std::vector<int>& row_blocks, - const std::vector<int>& col_blocks, - std::vector<int>* ordering); - - // Find a fill reducing approximate minimum degree - // ordering. ordering is expected to be large enough to hold the - // ordering. - bool ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering); - - // Find a fill reducing ordering using nested dissection. - bool NestedDissectionOrdering(cholmod_sparse* matrix, int* ordering); + // fill-reducing ordering algorithms can either be run on the block + // or the scalar form of these matrices. But since the underlying + // matrices are block oriented, it is worth running the fill + // reducing ordering on just the block structure of these matrices + // and then lifting these block orderings to a full scalar + // ordering. This preserves the block structure of the permuted + // matrix, and exposes more of the super-nodal structure of the + // matrix to the numerical factorization routines. + bool BlockOrdering(const cholmod_sparse* A, + OrderingType ordering_type, + const std::vector<int>& row_blocks, + const std::vector<int>& col_blocks, + std::vector<int>* ordering); // Nested dissection is only available if SuiteSparse is compiled // with Metis support. - static bool IsNestedDissectionAvailable() { -#ifdef CERES_NO_METIS - return false; -#else - return true; -#endif - } + static bool IsNestedDissectionAvailable(); // Find a fill reducing approximate minimum degree // ordering. constraints is an array which associates with each @@ -300,6 +283,9 @@ class CERES_NO_EXPORT SuiteSparse { public: + // Nested dissection is only available if SuiteSparse is compiled + // with Metis support. + static bool IsNestedDissectionAvailable() { return false; } void Free(void* /*arg*/) {} };
diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc index daf6601..9467b5a 100644 --- a/internal/ceres/trust_region_preprocessor.cc +++ b/internal/ceres/trust_region_preprocessor.cc
@@ -210,7 +210,10 @@ 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.use_postordering = options.use_postordering; + 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)) { @@ -236,7 +239,7 @@ // TODO(sameeragarwal): Implement the reordering of parameter // blocks for CX_SPARSE. if (options.sparse_linear_algebra_library_type == CX_SPARSE) { - pp->linear_solver_options.use_postordering = true; + pp->linear_solver_options.ordering_type = OrderingType::AMD; } } }
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc index f04442d..bde2814 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -106,14 +106,7 @@ LinearSolver::Options sparse_cholesky_options; sparse_cholesky_options.sparse_linear_algebra_library_type = options_.sparse_linear_algebra_library_type; - - // The preconditioner's sparsity is not available in the - // preprocessor, so the columns of the Jacobian have not been - // reordered to minimize fill in when computing its sparse Cholesky - // factorization. So we must tell the SparseCholesky object to - // perform approximate minimum-degree reordering, which is done by - // setting use_postordering to true. - sparse_cholesky_options.use_postordering = true; + sparse_cholesky_options.ordering_type = options.ordering_type; sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options); const time_t init_time = time(nullptr);