diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 8fff9f0bc..59420c91a 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -96,9 +96,6 @@ DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the " "accuracy of each linear solve of the truncated newton step. " "Changing this parameter can affect solve performance."); -DEFINE_bool(use_block_amd, true, "Use a block oriented fill reducing " - "ordering."); - DEFINE_int32(num_threads, 1, "Number of threads."); DEFINE_int32(num_iterations, 5, "Number of iterations."); DEFINE_double(max_solver_time, 1e32, "Maximum solve time in seconds."); @@ -139,8 +136,6 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) { const int camera_block_size = bal_problem->camera_block_size(); double* cameras = bal_problem->mutable_cameras(); - options->use_block_amd = FLAGS_use_block_amd; - if (options->use_inner_iterations) { if (FLAGS_blocks_for_inner_iterations == "cameras") { LOG(INFO) << "Camera blocks for inner iterations"; diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 8c2ff32d8..97d082d80 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -95,13 +95,8 @@ class Solver { #endif num_linear_solver_threads = 1; - -#if defined(CERES_NO_SUITESPARSE) - use_block_amd = false; -#else - use_block_amd = true; -#endif linear_solver_ordering = NULL; + use_postordering = false; use_inner_iterations = false; inner_iteration_ordering = NULL; linear_solver_min_num_iterations = 1; @@ -356,19 +351,29 @@ class Solver { // deallocate the memory when destroyed. ParameterBlockOrdering* linear_solver_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 (for - // SPARSE_NORMAL_CHOLESKY, SPARSE_SCHUR and ITERATIVE in - // conjunction with CLUSTER_TRIDIAGONAL AND CLUSTER_JACOBI - // preconditioners), the fill-reducing ordering algorithms can - // either be run on the block or the scalar form of these matrices. - // Running it on the block form exposes more of the super-nodal - // structure of the matrix to the factorization routines. Setting - // this parameter to true runs the ordering algorithms in block - // form. Currently this option only makes sense with - // sparse_linear_algebra_library = SUITE_SPARSE. - bool use_block_amd; + // Note: This option only applies to the SPARSE_NORMAL_CHOLESKY + // solver when used with SUITE_SPARSE. + + // Sparse Cholesky factorization algorithms use a fill-reducing + // ordering to permute the columns of the Jacobian matrix. There + // are two ways of doing this. + + // 1. Compute the Jacobian matrix in some order and then have the + // factorization algorithm permute the columns of the Jacobian. + + // 2. Compute the Jacobian with its columns already permuted. + + // The first option incurs a significant memory penalty. The + // factorization algorithm has to make a copy of the permuted + // Jacobian matrix, thus Ceres pre-permutes the columns of the + // Jacobian matrix and generally speaking, there is no performance + // penalty for doing so. + + // 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 + // // matrix. Setting use_postordering to true enables this tradeoff. + bool use_postordering; // Some non-linear least squares problems have additional // structure in the way the parameter blocks interact that it is diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc index cf5562ac7..15e0bdcd8 100644 --- a/internal/ceres/iterative_schur_complement_solver.cc +++ b/internal/ceres/iterative_schur_complement_solver.cc @@ -99,7 +99,6 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl( preconditioner_options.type = options_.preconditioner_type; preconditioner_options.sparse_linear_algebra_library = options_.sparse_linear_algebra_library; - preconditioner_options.use_block_amd = options_.use_block_amd; preconditioner_options.num_threads = options_.num_threads; preconditioner_options.row_block_size = options_.row_block_size; preconditioner_options.e_block_size = options_.e_block_size; diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index f4bd0fb6f..ca10faa24 100644 --- a/internal/ceres/linear_solver.h +++ b/internal/ceres/linear_solver.h @@ -74,7 +74,7 @@ class LinearSolver { : type(SPARSE_NORMAL_CHOLESKY), preconditioner_type(JACOBI), sparse_linear_algebra_library(SUITE_SPARSE), - use_block_amd(true), + use_postordering(false), min_num_iterations(1), max_num_iterations(1), num_threads(1), @@ -90,8 +90,8 @@ class LinearSolver { SparseLinearAlgebraLibraryType sparse_linear_algebra_library; - // See solver.h for explanation of this option. - bool use_block_amd; + // See solver.h for information about this flag. + bool use_postordering; // Number of internal iterations that the solver uses. This // parameter only makes sense for iterative solvers like CG. diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h index bfc8464db..d7c882936 100644 --- a/internal/ceres/preconditioner.h +++ b/internal/ceres/preconditioner.h @@ -47,7 +47,6 @@ class Preconditioner : public LinearOperator { Options() : type(JACOBI), sparse_linear_algebra_library(SUITE_SPARSE), - use_block_amd(true), num_threads(1), row_block_size(Eigen::Dynamic), e_block_size(Eigen::Dynamic), @@ -58,9 +57,6 @@ class Preconditioner : public LinearOperator { SparseLinearAlgebraLibraryType sparse_linear_algebra_library; - // See solver.h for explanation of this option. - bool use_block_amd; - // If possible, how many threads the preconditioner can use. int num_threads; diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index 875f45510..8afb12150 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc @@ -286,11 +286,7 @@ bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse( // Symbolic factorization is computed if we don't already have one handy. if (factor_ == NULL) { - if (options().use_block_amd) { - factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_); - } else { - factor_ = ss_.AnalyzeCholesky(cholmod_lhs); - } + factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_); } cholmod_dense* cholmod_solution = diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index 0d7b83c40..8217484e6 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc @@ -1004,7 +1004,8 @@ Program* SolverImpl::CreateReducedProgram(Solver::Options* options, return transformed_program.release(); } - if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY) { + if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY && + options->sparse_linear_algebra_library == SUITE_SPARSE) { ReorderProgramForSparseNormalCholesky(transformed_program.get()); return transformed_program.release(); } @@ -1091,11 +1092,10 @@ LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options, linear_solver_options.preconditioner_type = options->preconditioner_type; linear_solver_options.sparse_linear_algebra_library = options->sparse_linear_algebra_library; - + linear_solver_options.use_postordering = options->use_postordering; linear_solver_options.num_threads = options->num_linear_solver_threads; options->num_linear_solver_threads = linear_solver_options.num_threads; - linear_solver_options.use_block_amd = options->use_block_amd; const map >& groups = options->linear_solver_ordering->group_to_elements(); for (map >::const_iterator it = groups.begin(); diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc index 0cd064b55..bc1f98334 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/internal/ceres/sparse_normal_cholesky_solver.cc @@ -204,14 +204,15 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( event_logger.AddEvent("Setup"); if (factor_ == NULL) { - if (options_.use_block_amd) { + if (options_.use_postordering) { factor_ = ss_.BlockAnalyzeCholesky(&lhs, A->col_blocks(), A->row_blocks()); } else { - factor_ = ss_.AnalyzeCholesky(&lhs); + factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs); } } + event_logger.AddEvent("Analysis"); cholmod_dense* sol = ss_.SolveCholesky(&lhs, factor_, rhs); diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc index 0b0d593cf..c8a15c026 100644 --- a/internal/ceres/unsymmetric_linear_solver_test.cc +++ b/internal/ceres/unsymmetric_linear_solver_test.cc @@ -62,7 +62,6 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { LinearSolver::Options options; options.type = linear_solver_type; options.sparse_linear_algebra_library = sparse_linear_algebra_library; - options.use_block_amd = false; scoped_ptr solver(LinearSolver::Create(options)); LinearSolver::PerSolveOptions per_solve_options; diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc index f5c46bfcc..4b1e26af6 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc @@ -421,11 +421,7 @@ bool VisibilityBasedPreconditioner::Factorize() { // Symbolic factorization is computed if we don't already have one handy. if (factor_ == NULL) { - if (options_.use_block_amd) { - factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_); - } else { - factor_ = ss_.AnalyzeCholesky(lhs); - } + factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_); } bool status = ss_.Cholesky(lhs, factor_);