diff --git a/docs/source/solving.rst b/docs/source/solving.rst index f574ec321..f17c695cd 100644 --- a/docs/source/solving.rst +++ b/docs/source/solving.rst @@ -1853,6 +1853,7 @@ elimination group [LiSaad]_. TrustRegionStrategyType trust_region_strategy_type; DoglegType dogleg_type; + DenseLinearAlgebraLibraryType dense_linear_algebra_library_type; SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; LineSearchDirectionType line_search_direction_type; diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 358650135..224ad748e 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -132,7 +132,6 @@ void SetLinearSolver(Solver::Options* options) { FLAGS_dense_linear_algebra_library, &options->dense_linear_algebra_library_type)); options->num_linear_solver_threads = FLAGS_num_threads; - } void SetOrdering(BALProblem* bal_problem, Solver::Options* options) { diff --git a/include/ceres/solver.h b/include/ceres/solver.h index ab1d35009..25b762a7b 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -99,12 +99,13 @@ class Solver { preconditioner_type = JACOBI; + dense_linear_algebra_library_type = EIGEN; sparse_linear_algebra_library_type = SUITE_SPARSE; #if defined(CERES_NO_SUITESPARSE) && !defined(CERES_NO_CXSPARSE) sparse_linear_algebra_library_type = CX_SPARSE; #endif - dense_linear_algebra_library_type = EIGEN; + num_linear_solver_threads = 1; linear_solver_ordering = NULL; use_postordering = false; @@ -384,12 +385,6 @@ class Solver { // Type of preconditioner to use with the iterative linear solvers. PreconditionerType preconditioner_type; - // Ceres supports using multiple sparse linear algebra libraries - // for sparse matrix ordering and factorizations. Currently, - // SUITE_SPARSE and CX_SPARSE are the valid choices, depending on - // whether they are linked into Ceres at build time. - SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; - // Ceres supports using multiple dense linear algebra libraries // for dense matrix factorizations. Currently EIGEN and LAPACK are // the valid choices. EIGEN is always available, LAPACK refers to @@ -403,6 +398,12 @@ class Solver { // performance. DenseLinearAlgebraLibraryType dense_linear_algebra_library_type; + // Ceres supports using multiple sparse linear algebra libraries + // for sparse matrix ordering and factorizations. Currently, + // SUITE_SPARSE and CX_SPARSE are the valid choices, depending on + // whether they are linked into Ceres at build time. + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; + // Number of threads used by Ceres to solve the Newton // step. Currently only the SPARSE_SCHUR solver is capable of // using this setting. diff --git a/internal/ceres/blas.cc b/internal/ceres/blas.cc index d63ca5159..f79b1ebfa 100644 --- a/internal/ceres/blas.cc +++ b/internal/ceres/blas.cc @@ -72,7 +72,7 @@ void BLAS::SymmetricRankKUpdate(int num_rows, c, &ldc); #endif -}; +} } // namespace internal } // namespace ceres diff --git a/internal/ceres/blas.h b/internal/ceres/blas.h index a5c6862f1..2ab666395 100644 --- a/internal/ceres/blas.h +++ b/internal/ceres/blas.h @@ -38,7 +38,6 @@ namespace internal { class BLAS { public: - // transpose = true : c = alpha * a'a + beta * c; // transpose = false : c = alpha * aa' + beta * c; // diff --git a/internal/ceres/compressed_col_sparse_matrix_utils.h b/internal/ceres/compressed_col_sparse_matrix_utils.h index 494a8a5c8..c8de2a159 100644 --- a/internal/ceres/compressed_col_sparse_matrix_utils.h +++ b/internal/ceres/compressed_col_sparse_matrix_utils.h @@ -80,7 +80,7 @@ void SolveUpperTriangularInPlace(IntegerType num_cols, rhs_and_solution[r] -= v * rhs_and_solution[c]; } } -}; +} // Solve the linear system // @@ -100,8 +100,8 @@ void SolveUpperTriangularTransposeInPlace(IntegerType num_cols, rhs_and_solution[c] -= v * rhs_and_solution[r]; } rhs_and_solution[c] = rhs_and_solution[c] / values[cols[c + 1] - 1]; - }; -}; + } +} // Given a upper triangular matrix R in compressed column form, solve // the linear system, @@ -131,10 +131,10 @@ void SolveRTRWithSparseRHS(IntegerType num_cols, solution[c] -= v * solution[r]; } solution[c] = solution[c] / values[cols[c + 1] - 1]; - }; - SolveUpperTriangularInPlace(num_cols, rows, cols, values, solution); -}; + } + SolveUpperTriangularInPlace(num_cols, rows, cols, values, solution); +} } // namespace internal } // namespace ceres diff --git a/internal/ceres/dense_normal_cholesky_solver.cc b/internal/ceres/dense_normal_cholesky_solver.cc index 677987e05..fbf3cbec9 100644 --- a/internal/ceres/dense_normal_cholesky_solver.cc +++ b/internal/ceres/dense_normal_cholesky_solver.cc @@ -54,11 +54,11 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveImpl( const double* b, const LinearSolver::PerSolveOptions& per_solve_options, double* x) { - if (options_.dense_linear_algebra_library_type == EIGEN) { - return SolveUsingEigen(A, b, per_solve_options, x); - } else { - return SolveUsingLAPACK(A, b, per_solve_options, x); - } + if (options_.dense_linear_algebra_library_type == EIGEN) { + return SolveUsingEigen(A, b, per_solve_options, x); + } else { + return SolveUsingLAPACK(A, b, per_solve_options, x); + } } LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen( @@ -93,7 +93,6 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen( } event_logger.AddEvent("Product"); - // Use dsyrk instead for the product. LinearSolver::Summary summary; summary.num_iterations = 1; summary.termination_type = TOLERANCE; @@ -139,7 +138,8 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingLAPACK( // TODO(sameeragarwal): Replace this with a gemv call for true blasness. // rhs = A'b - VectorRef(x, num_cols) = A->matrix().transpose() * ConstVectorRef(b, A->num_rows()); + VectorRef(x, num_cols) = + A->matrix().transpose() * ConstVectorRef(b, A->num_rows()); event_logger.AddEvent("Product"); const int info = LAPACK::SolveInPlaceUsingCholesky(num_cols, lhs.data(), x); diff --git a/internal/ceres/dense_qr_solver.cc b/internal/ceres/dense_qr_solver.cc index c61a96fa5..d76d58b51 100644 --- a/internal/ceres/dense_qr_solver.cc +++ b/internal/ceres/dense_qr_solver.cc @@ -76,6 +76,8 @@ LinearSolver::Summary DenseQRSolver::SolveUsingLAPACK( A->AppendDiagonal(per_solve_options.D); } + // TODO(sameeragarwal): Since we are copying anyways, the diagonal + // can be appended to the matrix instead of doing it on A. lhs_ = A->matrix(); if (per_solve_options.D != NULL) { @@ -91,7 +93,8 @@ LinearSolver::Summary DenseQRSolver::SolveUsingLAPACK( rhs_.head(num_rows) = ConstVectorRef(b, num_rows); if (work_.rows() == 1) { - const int work_size = LAPACK::EstimateWorkSizeForQR(lhs_.rows(), lhs_.cols()); + const int work_size = + LAPACK::EstimateWorkSizeForQR(lhs_.rows(), lhs_.cols()); VLOG(3) << "Working memory for Dense QR factorization: " << work_size * sizeof(double); work_.resize(work_size); @@ -128,7 +131,6 @@ LinearSolver::Summary DenseQRSolver::SolveUsingEigen( const int num_rows = A->num_rows(); const int num_cols = A->num_cols(); - if (per_solve_options.D != NULL) { // Temporarily append a diagonal block to the A matrix, but undo // it before returning the matrix to the user. diff --git a/internal/ceres/lapack.cc b/internal/ceres/lapack.cc index 6e7d21b58..73bfa69cb 100644 --- a/internal/ceres/lapack.cc +++ b/internal/ceres/lapack.cc @@ -101,11 +101,22 @@ int LAPACK::EstimateWorkSizeForQR(int num_rows, int num_cols) { int lwork = -1; double work; int info = 0; - dgels_(&trans, &num_rows, &num_cols, &nrhs, NULL, &num_rows, NULL, &num_rows, &work, &lwork, &info); + dgels_(&trans, + &num_rows, + &num_cols, + &nrhs, + NULL, + &num_rows, + NULL, + &num_rows, + &work, + &lwork, + &info); + CHECK_EQ(info, 0); return work; #endif -}; +} int LAPACK::SolveUsingQR(int num_rows, int num_cols, @@ -126,10 +137,21 @@ int LAPACK::SolveUsingQR(int num_rows, int info = 0; double* lhs = const_cast(in_lhs); - dgels_(&trans, &m, &n, &nrhs, lhs, &lda, rhs_and_solution, &ldb, work, &work_size, &info); + dgels_(&trans, + &m, + &n, + &nrhs, + lhs, + &lda, + rhs_and_solution, + &ldb, + work, + &work_size, + &info); + return info; #endif -}; +} } // namespace internal } // namespace ceres diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc index 745ea8ec3..d91c16228 100644 --- a/internal/ceres/schur_complement_solver_test.cc +++ b/internal/ceres/schur_complement_solver_test.cc @@ -160,7 +160,8 @@ TEST_F(SchurComplementSolverTest, LAPACKBasedDenseSchurWithLargeProblem) { #ifndef CERES_NO_SUITESPARSE TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) { - ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); + ComputeAndCompareSolutions( + 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); } @@ -172,7 +173,8 @@ TEST_F(SchurComplementSolverTest, TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) { - ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); + ComputeAndCompareSolutions( + 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false); } diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index bd731feb4..06c65f921 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc @@ -1231,7 +1231,7 @@ LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options, // done. #if !defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CAMD) if (IsSchurType(linear_solver_options.type) && - linear_solver_options.sparse_linear_algebra_library_type == SUITE_SPARSE) { + options->sparse_linear_algebra_library_type == SUITE_SPARSE) { linear_solver_options.use_postordering = true; } #endif diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h index c1772924d..859e20e41 100644 --- a/internal/ceres/suitesparse.h +++ b/internal/ceres/suitesparse.h @@ -268,7 +268,7 @@ class SuiteSparse { } // namespace internal } // namespace ceres -#else // CERES_NO_SUITESPARSE +#else // CERES_NO_SUITESPARSE class SuiteSparse {}; typedef void cholmod_factor;