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);