diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index d03c2ff59..87e18e110 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -95,7 +95,7 @@ DEFINE_string(sparse_linear_algebra_library, "suite_sparse", "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_int32(max_num_refinement_iterations, 0, "Iterative refinement iterations" 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 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) { // 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 026fc1c08..43146b110 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -395,12 +395,33 @@ class CERES_EXPORT Solver { #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 @@ class CERES_EXPORT Solver { // 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 @@ class CERES_EXPORT Solver { // {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 @@ class CERES_EXPORT Solver { // 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 @@ class CERES_EXPORT Solver { 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 linear_solver_ordering_given; diff --git a/include/ceres/types.h b/include/ceres/types.h index e52242381..2fc5abccd 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -183,6 +183,30 @@ enum SparseLinearAlgebraLibraryType { 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 const char* SparseLinearAlgebraLibraryTypeToString( 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 c6765fef7..1f50a39ef 100644 --- a/internal/ceres/reorder_program.cc +++ b/internal/ceres/reorder_program.cc @@ -85,7 +85,6 @@ static int MinParameterBlock(const ResidualBlock* residual_block, return min_parameter_block_position; } -#if defined(CERES_USE_EIGEN_SPARSE) Eigen::SparseMatrix CreateBlockJacobian( const TripletSparseMatrix& block_jacobian_transpose) { using SparseMatrix = Eigen::SparseMatrix; @@ -105,9 +104,9 @@ Eigen::SparseMatrix CreateBlockJacobian( 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& parameter_blocks, const ParameterBlockOrdering& parameter_block_ordering, @@ -120,40 +119,55 @@ void OrderingForSparseNormalCholeskyUsingSuiteSparse( cholmod_sparse* block_jacobian_transpose = ss.CreateSparseMatrix( const_cast(&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 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 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 OrderingForSparseNormalCholeskyUsingCXSparse( } 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 @@ bool LexicographicallyOrderResidualBlocks( // 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 @@ static void MaybeReorderSchurComplementColumnsUsingSuiteSparse( #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 @@ static void MaybeReorderSchurComplementColumnsUsingEigen( 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 @@ bool ReorderProgramForSchurTypeLinearSolver( 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 ReorderProgramForSchurTypeLinearSolver( 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 @@ bool ReorderProgramForSparseCholesky( 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 @@ bool ReorderProgramForSparseCholesky( } 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 @@ int ReorderResidualBlocksByPartition( 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 64f409504..8ae456ec4 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 LexicographicallyOrderResidualBlocks( 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 @@ CERES_NO_EXPORT bool ReorderProgramForSchurTypeLinearSolver( // 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 @@ CERES_NO_EXPORT int ReorderResidualBlocksByPartition( const std::unordered_set& 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 0fbab4e10..872f3812c 100644 --- a/internal/ceres/reorder_program_test.cc +++ b/internal/ceres/reorder_program_test.cc @@ -185,6 +185,7 @@ class ReorderProgramForSparseCholeskyUsingSuiteSparseTest 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 150c5550f..dbaba3e89 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 @@ bool CommonOptionsAreValid(const Solver::Options& options, string* error) { 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 @@ bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) { 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 @@ bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) { 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 @@ void PreSolveSummarize(const Solver::Options& options, &(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 @@ void PreSolveSummarize(const Solver::Options& options, 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 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, 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 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, internal::SetSummaryFinalCost(summary); - if (pp.reduced_program.get() != nullptr) { + if (pp.reduced_program != nullptr) { SummarizeReducedProgram(*pp.reduced_program, summary); } @@ -389,7 +417,7 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, // 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& evaluator_statistics = pp.evaluator->Statistics(); { @@ -411,7 +439,7 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, // 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& linear_solver_statistics = pp.linear_solver->Statistics(); const CallStatistics& call_stats = FindWithDefault( @@ -666,17 +694,6 @@ string Solver::Summary::FullReport() const { 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 @@ string Solver::Summary::FullReport() const { 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 9467b5ac6..ed62aa9d8 100644 --- a/internal/ceres/trust_region_preprocessor.cc +++ b/internal/ceres/trust_region_preprocessor.cc @@ -113,6 +113,7 @@ bool ReorderProgram(PreprocessedProblem* pp) { 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 @@ bool ReorderProgram(PreprocessedProblem* pp) { !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 @@ bool ReorderProgram(PreprocessedProblem* pp) { 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 @@ bool SetupLinearSolver(PreprocessedProblem* pp) { 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 @@ bool SetupLinearSolver(PreprocessedProblem* pp) { 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 @@ bool SetupLinearSolver(PreprocessedProblem* pp) { 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 48242678b..eea031237 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -124,6 +124,23 @@ bool StringToSparseLinearAlgebraLibraryType( 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 bde281406..0a7b19f97 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc @@ -106,7 +106,7 @@ VisibilityBasedPreconditioner::VisibilityBasedPreconditioner( 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);