diff options
Diffstat (limited to 'extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc')
| -rw-r--r-- | extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc | 250 |
1 files changed, 207 insertions, 43 deletions
diff --git a/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc b/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc index cf5bb235b46..94f7e5803c4 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc @@ -28,11 +28,6 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARSE) - #include "ceres/sparse_normal_cholesky_solver.h" #include <algorithm> @@ -48,9 +43,67 @@ #include "ceres/triplet_sparse_matrix.h" #include "ceres/types.h" #include "ceres/wall_time.h" +#include "Eigen/SparseCore" + +#ifdef CERES_USE_EIGEN_SPARSE +#include "Eigen/SparseCholesky" +#endif namespace ceres { namespace internal { +namespace { + +#ifdef CERES_USE_EIGEN_SPARSE +// A templated factorized and solve function, which allows us to use +// the same code independent of whether a AMD or a Natural ordering is +// used. +template <typename SimplicialCholeskySolver> +LinearSolver::Summary SimplicialLDLTSolve( + Eigen::MappedSparseMatrix<double, Eigen::ColMajor>& lhs, + const bool do_symbolic_analysis, + SimplicialCholeskySolver* solver, + double* rhs_and_solution, + EventLogger* event_logger) { + LinearSolver::Summary summary; + summary.num_iterations = 1; + summary.termination_type = LINEAR_SOLVER_SUCCESS; + summary.message = "Success."; + + if (do_symbolic_analysis) { + solver->analyzePattern(lhs); + event_logger->AddEvent("Analyze"); + if (solver->info() != Eigen::Success) { + summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; + summary.message = + "Eigen failure. Unable to find symbolic factorization."; + return summary; + } + } + + solver->factorize(lhs); + event_logger->AddEvent("Factorize"); + if (solver->info() != Eigen::Success) { + summary.termination_type = LINEAR_SOLVER_FAILURE; + summary.message = "Eigen failure. Unable to find numeric factorization."; + return summary; + } + + const Vector rhs = VectorRef(rhs_and_solution, lhs.cols()); + + VectorRef(rhs_and_solution, lhs.cols()) = solver->solve(rhs); + event_logger->AddEvent("Solve"); + if (solver->info() != Eigen::Success) { + summary.termination_type = LINEAR_SOLVER_FAILURE; + summary.message = "Eigen failure. Unable to do triangular solve."; + return summary; + } + + return summary; +} + +#endif // CERES_USE_EIGEN_SPARSE + +} // namespace SparseNormalCholeskySolver::SparseNormalCholeskySolver( const LinearSolver::Options& options) @@ -60,19 +113,15 @@ SparseNormalCholeskySolver::SparseNormalCholeskySolver( } void SparseNormalCholeskySolver::FreeFactorization() { -#ifndef CERES_NO_SUITESPARSE if (factor_ != NULL) { ss_.Free(factor_); factor_ = NULL; } -#endif // CERES_NO_SUITESPARSE -#ifndef CERES_NO_CXSPARSE if (cxsparse_factor_ != NULL) { cxsparse_.Free(cxsparse_factor_); cxsparse_factor_ = NULL; } -#endif // CERES_NO_CXSPARSE } SparseNormalCholeskySolver::~SparseNormalCholeskySolver() { @@ -111,6 +160,9 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl( case CX_SPARSE: summary = SolveImplUsingCXSparse(A, per_solve_options, x); break; + case EIGEN_SPARSE: + summary = SolveImplUsingEigen(A, per_solve_options, x); + break; default: LOG(FATAL) << "Unknown sparse linear algebra library : " << options_.sparse_linear_algebra_library_type; @@ -123,11 +175,120 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl( return summary; } -#ifndef CERES_NO_CXSPARSE +LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingEigen( + CompressedRowSparseMatrix* A, + const LinearSolver::PerSolveOptions& per_solve_options, + double * rhs_and_solution) { +#ifndef CERES_USE_EIGEN_SPARSE + + LinearSolver::Summary summary; + summary.num_iterations = 0; + summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; + summary.message = + "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."; + return summary; + +#else + + EventLogger event_logger("SparseNormalCholeskySolver::Eigen::Solve"); + // Compute the normal equations. J'J delta = J'f and solve them + // using a sparse Cholesky factorization. Notice that when compared + // to SuiteSparse we have to explicitly compute the normal equations + // before they can be factorized. CHOLMOD/SuiteSparse on the other + // hand can just work off of Jt to compute the Cholesky + // factorization of the normal equations. + // + // TODO(sameeragarwal): See note about how this maybe a bad idea for + // dynamic sparsity. + if (outer_product_.get() == NULL) { + outer_product_.reset( + CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( + *A, &pattern_)); + } + + CompressedRowSparseMatrix::ComputeOuterProduct( + *A, pattern_, outer_product_.get()); + + // Map to an upper triangular column major matrix. + // + // outer_product_ is a compressed row sparse matrix and in lower + // triangular form, when mapped to a compressed column sparse + // matrix, it becomes an upper triangular matrix. + Eigen::MappedSparseMatrix<double, Eigen::ColMajor> AtA( + outer_product_->num_rows(), + outer_product_->num_rows(), + outer_product_->num_nonzeros(), + outer_product_->mutable_rows(), + outer_product_->mutable_cols(), + outer_product_->mutable_values()); + + // Dynamic sparsity means that we cannot depend on a static analysis + // of sparsity structure of the jacobian, so we compute a new + // symbolic factorization every time. + if (options_.dynamic_sparsity) { + amd_ldlt_.reset(NULL); + } + + bool do_symbolic_analysis = false; + + // If using post ordering or dynamic sparsity, or an old version of + // Eigen, we cannot depend on a preordered jacobian, so we work with + // a SimplicialLDLT decomposition with AMD ordering. + if (options_.use_postordering || + options_.dynamic_sparsity || + !EIGEN_VERSION_AT_LEAST(3, 2, 2)) { + if (amd_ldlt_.get() == NULL) { + amd_ldlt_.reset(new SimplicialLDLTWithAMDOrdering); + do_symbolic_analysis = true; + } + + return SimplicialLDLTSolve(AtA, + do_symbolic_analysis, + amd_ldlt_.get(), + rhs_and_solution, + &event_logger); + } + +#if EIGEN_VERSION_AT_LEAST(3,2,2) + // The common case + if (natural_ldlt_.get() == NULL) { + natural_ldlt_.reset(new SimplicialLDLTWithNaturalOrdering); + do_symbolic_analysis = true; + } + + return SimplicialLDLTSolve(AtA, + do_symbolic_analysis, + natural_ldlt_.get(), + rhs_and_solution, + &event_logger); +#endif + +#endif // EIGEN_USE_EIGEN_SPARSE +} + + + LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( CompressedRowSparseMatrix* A, const LinearSolver::PerSolveOptions& per_solve_options, double * rhs_and_solution) { +#ifdef CERES_NO_CXSPARSE + + LinearSolver::Summary summary; + summary.num_iterations = 0; + summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; + summary.message = + "SPARSE_NORMAL_CHOLESKY cannot be used with CX_SPARSE " + "because Ceres was not built with support for CXSparse. " + "This requires enabling building with -DCXSPARSE=ON."; + + return summary; + +#else + EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve"); LinearSolver::Summary summary; @@ -137,11 +298,14 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( // Compute the normal equations. J'J delta = J'f and solve them // using a sparse Cholesky factorization. Notice that when compared - // to SuiteSparse we have to explicitly compute the transpose of Jt, - // and then the normal equations before they can be - // factorized. CHOLMOD/SuiteSparse on the other hand can just work - // off of Jt to compute the Cholesky factorization of the normal - // equations. + // to SuiteSparse we have to explicitly compute the normal equations + // before they can be factorized. CHOLMOD/SuiteSparse on the other + // hand can just work off of Jt to compute the Cholesky + // factorization of the normal equations. + // + // TODO(sameeragarwal): If dynamic sparsity is enabled, then this is + // not a good idea performance wise, since the jacobian has far too + // many entries and the program will go crazy with memory. if (outer_product_.get() == NULL) { outer_product_.reset( CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( @@ -179,30 +343,35 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; summary.message = "CXSparse failure. Unable to find symbolic factorization."; - } else if (!cxsparse_.SolveCholesky(AtA, cxsparse_factor_, rhs_and_solution)) { + } else if (!cxsparse_.SolveCholesky(AtA, + cxsparse_factor_, + rhs_and_solution)) { summary.termination_type = LINEAR_SOLVER_FAILURE; + summary.message = "CXSparse::SolveCholesky failed."; } event_logger.AddEvent("Solve"); return summary; -} -#else -LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( - CompressedRowSparseMatrix* A, - const LinearSolver::PerSolveOptions& per_solve_options, - double * rhs_and_solution) { - LOG(FATAL) << "No CXSparse support in Ceres."; - - // Unreachable but MSVC does not know this. - return LinearSolver::Summary(); -} #endif +} -#ifndef CERES_NO_SUITESPARSE LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, const LinearSolver::PerSolveOptions& per_solve_options, double * rhs_and_solution) { +#ifdef CERES_NO_SUITESPARSE + + LinearSolver::Summary summary; + summary.num_iterations = 0; + summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; + summary.message = + "SPARSE_NORMAL_CHOLESKY cannot be used with SUITE_SPARSE " + "because Ceres was not built with support for SuiteSparse. " + "This requires enabling building with -DSUITESPARSE=ON."; + return summary; + +#else + EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve"); LinearSolver::Summary summary; summary.termination_type = LINEAR_SOLVER_SUCCESS; @@ -226,7 +395,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( if (options_.dynamic_sparsity) { factor_ = ss_.AnalyzeCholesky(&lhs, &summary.message); } else { - factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message); + factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, + &summary.message); } } } @@ -234,6 +404,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( if (factor_ == NULL) { summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; + // No need to set message as it has already been set by the + // symbolic analysis routines above. return summary; } @@ -242,7 +414,9 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( return summary; } - cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols); + cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution, + num_cols, + num_cols); cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message); event_logger.AddEvent("Solve"); @@ -251,25 +425,15 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( memcpy(rhs_and_solution, solution->x, num_cols * sizeof(*rhs_and_solution)); ss_.Free(solution); } else { + // No need to set message as it has already been set by the + // numeric factorization routine above. summary.termination_type = LINEAR_SOLVER_FAILURE; } event_logger.AddEvent("Teardown"); return summary; -} -#else -LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( - CompressedRowSparseMatrix* A, - const LinearSolver::PerSolveOptions& per_solve_options, - double * rhs_and_solution) { - LOG(FATAL) << "No SuiteSparse support in Ceres."; - - // Unreachable but MSVC does not know this. - return LinearSolver::Summary(); -} #endif +} } // namespace internal } // namespace ceres - -#endif // !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARSE) |