Multiple sparse linear algebra backends. 1. Added support for CXSparse - SparseNormalCholesky and SchurComplementSolver support SuiteSparse and CXSparse now. I am not sure I will add suport for visibility based preconditioning using CXSparse. Its not a high priority. 2. New enum SparseLinearAlgebraLibraryType which allows the user to indicate which sparse linear algebra library should be used. 3. Updated tests for SolverImpl and system_test. 4. Build system changes to automatically detect CXSparse and link to it by default -- just like SuiteSparse. 5. Minor bug fixes dealing in the cmake files and VBP. 6. Changed the order of the system test. 7. Deduped the unsymmetric linear solver test. Change-Id: I33252a103c87b722ecb7ed7b5f0ae7fd91249244
diff --git a/CMakeLists.txt b/CMakeLists.txt index 39d4c64..ab00deb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt
@@ -59,27 +59,31 @@ # Locations to search for Eigen SET(EIGEN_SEARCH_HEADERS ${SEARCH_HEADERS}) LIST(APPEND EIGEN_SEARCH_HEADERS /usr/include/eigen3) # Ubuntu 10.04's default location. -LIST(APPEND EIGEN_SEARCH_HEADERS /usr/local/include/eigen3) +LIST(APPEND EIGEN_SEARCH_HEADERS /usr/local/include/eigen3) LIST(APPEND EIGEN_SEARCH_HEADERS /usr/local/homebrew/include/eigen3) # Mac OS X LIST(APPEND EIGEN_SEARCH_HEADERS /opt/local/var/macports/software/eigen3/opt/local/include/eigen3) # Mac OS X # Locations to search for SuiteSparse SET(SUITESPARSE_SEARCH_LIBS ${SEARCH_LIBS}) LIST(APPEND SUITESPARSE_SEARCH_LIBS /usr/lib/suitesparse) # Ubuntu -LIST(APPEND SUITESPARSE_SEARCH_LIBS /usr/local/lib/suitesparse) +LIST(APPEND SUITESPARSE_SEARCH_LIBS /usr/local/lib/suitesparse) LIST(APPEND SUITESPARSE_SEARCH_LIBS /opt/local/lib/ufsparse) # Mac OS X SET(SUITESPARSE_SEARCH_HEADERS ${SEARCH_HEADERS}) LIST(APPEND SUITESPARSE_SEARCH_HEADERS /usr/include/suitesparse) # Ubuntu -LIST(APPEND SUITESPARSE_SEARCH_HEADERS /usr/local/include/suitesparse) +LIST(APPEND SUITESPARSE_SEARCH_HEADERS /usr/local/include/suitesparse) LIST(APPEND SUITESPARSE_SEARCH_HEADERS /opt/local/include/ufsparse) # Mac OS X +SET(CXSPARSE_SEARCH_LIBS ${SEARCH_LIBS}) +SET(CXSPARSE_SEARCH_HEADERS ${SEARCH_HEADERS}) +LIST(APPEND CXSPARSE_SEARCH_HEADERS /usr/include/suitesparse) # Ubuntu + # Check for SuiteSparse dependencies MESSAGE("-- Check for AMD") SET(AMD_FOUND TRUE) FIND_LIBRARY(AMD_LIB NAMES amd PATHS ${SUITESPARSE_SEARCH_LIBS}) -IF (EXISTS ${AMD_LIB}) +IF (EXISTS ${AMD_LIB}) MESSAGE("-- Found AMD library: ${AMD_LIB}") ELSE (EXISTS ${AMD_LIB}) MESSAGE("-- Did not find AMD library") @@ -87,7 +91,7 @@ ENDIF (EXISTS ${AMD_LIB}) FIND_PATH(AMD_INCLUDE NAMES amd.h PATHS ${SUITESPARSE_SEARCH_HEADERS}) -IF (EXISTS ${AMD_INCLUDE}) +IF (EXISTS ${AMD_INCLUDE}) MESSAGE("-- Found AMD header in: ${AMD_INCLUDE}") ELSE (EXISTS ${AMD_INCLUDE}) MESSAGE("-- Did not find AMD header") @@ -98,7 +102,7 @@ SET(CAMD_FOUND TRUE) FIND_LIBRARY(CAMD_LIB NAMES camd PATHS ${SUITESPARSE_SEARCH_LIBS}) -IF (EXISTS ${CAMD_LIB}) +IF (EXISTS ${CAMD_LIB}) MESSAGE("-- Found CAMD library: ${CAMD_LIB}") ELSE (EXISTS ${CAMD_LIB}) MESSAGE("-- Did not find CAMD library") @@ -106,7 +110,7 @@ ENDIF (EXISTS ${CAMD_LIB}) FIND_PATH(CAMD_INCLUDE NAMES camd.h PATHS ${SUITESPARSE_SEARCH_HEADERS}) -IF (EXISTS ${CAMD_INCLUDE}) +IF (EXISTS ${CAMD_INCLUDE}) MESSAGE("-- Found CAMD header in: ${CAMD_INCLUDE}") ELSE (EXISTS ${CAMD_INCLUDE}) MESSAGE("-- Did not find CAMD header") @@ -117,7 +121,7 @@ SET(COLAMD_FOUND TRUE) FIND_LIBRARY(COLAMD_LIB NAMES colamd PATHS ${SUITESPARSE_SEARCH_LIBS}) -IF (EXISTS ${COLAMD_LIB}) +IF (EXISTS ${COLAMD_LIB}) MESSAGE("-- Found COLAMD library: ${COLAMD_LIB}") ELSE (EXISTS ${COLAMD_LIB}) MESSAGE("-- Did not find COLAMD library") @@ -125,7 +129,7 @@ ENDIF (EXISTS ${COLAMD_LIB}) FIND_PATH(COLAMD_INCLUDE NAMES colamd.h PATHS ${SUITESPARSE_SEARCH_HEADERS}) -IF (EXISTS ${COLAMD_INCLUDE}) +IF (EXISTS ${COLAMD_INCLUDE}) MESSAGE("-- Found COLAMD header in: ${COLAMD_INCLUDE}") ELSE (EXISTS ${COLAMD_INCLUDE}) MESSAGE("-- Did not find COLAMD header") @@ -136,7 +140,7 @@ SET(CCOLAMD_FOUND TRUE) FIND_LIBRARY(CCOLAMD_LIB NAMES ccolamd PATHS ${SUITESPARSE_SEARCH_LIBS}) -IF (EXISTS ${CCOLAMD_LIB}) +IF (EXISTS ${CCOLAMD_LIB}) MESSAGE("-- Found CCOLAMD library: ${CCOLAMD_LIB}") ELSE (EXISTS ${CCOLAMD_LIB}) MESSAGE("-- Did not find CCOLAMD library") @@ -144,7 +148,7 @@ ENDIF (EXISTS ${CCOLAMD_LIB}) FIND_PATH(CCOLAMD_INCLUDE NAMES ccolamd.h PATHS ${SUITESPARSE_SEARCH_HEADERS}) -IF (EXISTS ${CCOLAMD_INCLUDE}) +IF (EXISTS ${CCOLAMD_INCLUDE}) MESSAGE("-- Found CCOLAMD header in: ${CCOLAMD_INCLUDE}") ELSE (EXISTS ${CCOLAMD_INCLUDE}) MESSAGE("-- Did not find CCOLAMD header") @@ -155,7 +159,7 @@ SET(CHOLMOD_FOUND TRUE) FIND_LIBRARY(CHOLMOD_LIB NAMES cholmod PATHS ${SUITESPARSE_SEARCH_LIBS}) -IF (EXISTS ${CHOLMOD_LIB}) +IF (EXISTS ${CHOLMOD_LIB}) MESSAGE("-- Found CHOLMOD library: ${CHOLMOD_LIB}") ELSE (EXISTS ${CHOLMOD_LIB}) MESSAGE("-- Did not find CHOLMOD library") @@ -163,7 +167,7 @@ ENDIF (EXISTS ${CHOLMOD_LIB}) FIND_PATH(CHOLMOD_INCLUDE NAMES cholmod.h PATHS ${SUITESPARSE_SEARCH_HEADERS}) -IF (EXISTS ${CHOLMOD_INCLUDE}) +IF (EXISTS ${CHOLMOD_INCLUDE}) MESSAGE("-- Found CHOLMOD header in: ${CHOLMOD_INCLUDE}") ELSE (EXISTS ${CHOLMOD_INCLUDE}) MESSAGE("-- Did not find CHOLMOD header") @@ -173,7 +177,7 @@ MESSAGE("-- Check for METIS (optional)") FIND_LIBRARY(METIS_LIB NAMES metis PATHS ${SUITESPARSE_SEARCH_LIBS}) -IF (EXISTS ${METIS_LIB}) +IF (EXISTS ${METIS_LIB}) MESSAGE("-- Found METIS library: ${METIS_LIB}") ELSE (EXISTS ${METIS_LIB}) MESSAGE("-- Did not find METIS library") @@ -187,7 +191,7 @@ FIND_LIBRARY(LAPACK_LIB NAMES vecLib) ELSE (APPLE) FIND_LIBRARY(BLAS_LIB NAMES blas) - IF (EXISTS ${BLAS_LIB}) + IF (EXISTS ${BLAS_LIB}) MESSAGE("-- Found BLAS library: ${BLAS_LIB}") ELSE (EXISTS ${BLAS_LIB}) MESSAGE("-- Did not find BLAS library") @@ -196,14 +200,14 @@ FIND_LIBRARY(LAPACK_LIB NAMES lapack) ENDIF (APPLE) -IF (EXISTS ${LAPACK_LIB}) +IF (EXISTS ${LAPACK_LIB}) MESSAGE("-- Found LAPACK library: ${LAPACK_LIB}") ELSE (EXISTS ${LAPACK_LIB}) SET(BLAS_AND_LAPACK_FOUND FALSE) MESSAGE("-- Did not find LAPACK library") ENDIF (EXISTS ${LAPACK_LIB}) -SET(SUITESPARSE_FOUND +SET(SUITESPARSE_FOUND AMD_FOUND AND CAMD_FOUND AND COLAMD_FOUND AND @@ -215,20 +219,65 @@ # built with SuiteSparse support. -DSUITESPARSE=ON/OFF can be used to # enable/disable SuiteSparse explicitly. IF (DEFINED SUITESPARSE) - IF (SUITESPARSE AND NOT SUITESPARSE_FOUND) - MESSAGE(FATAL_ERROR "One or more of SuiteSparse's dependencies was not found") - ENDIF (SUITESPARSE AND NOT SUITESPARSE_FOUND) + IF (SUITESPARSE) + IF (NOT SUITESPARSE_FOUND) + MESSAGE(FATAL_ERROR "One or more of SuiteSparse's dependencies was not found") + ENDIF (NOT SUITESPARSE_FOUND) + ELSE (SUITESPARSE) + ADD_DEFINITIONS(-DCERES_NO_SUITESPARSE) + ENDIF (SUITESPARSE) ELSE (DEFINED SUITESPARSE) IF (SUITESPARSE_FOUND) MESSAGE("-- Found all SuiteSparse dependencies. Building with SuiteSparse") SET(SUITESPARSE ON) ELSE (SUITESPARSE_FOUND) - MESSAGE("-- Did not find all SuiteSparse dependencies. Building without SuiteSparse") + MESSAGE("-- Did not find all SuiteSparse dependencies. Building without SuiteSparse") SET(SUITESPARSE OFF) ADD_DEFINITIONS(-DCERES_NO_SUITESPARSE) ENDIF (SUITESPARSE_FOUND) ENDIF (DEFINED SUITESPARSE) +# By default, if all of CXSparse's dependencies are found, Ceres is +# built with CXSparse support. -DCXSPARSE=ON/OFF can be used to +# enable/disable CXSparse explicitly. +MESSAGE("-- Check for CXSparse") +SET(CXSPARSE_FOUND ON) + +FIND_LIBRARY(CXSPARSE_LIB NAMES cxsparse PATHS ${CXSPARSE_SEARCH_LIBS}) +IF (EXISTS ${CXSPARSE_LIB}) + MESSAGE("-- Found CXSparse library in: ${CXSPARSE_LIB}") +ELSE (EXISTS ${CXSPARSE_LIB}) + MESSAGE("-- Did not find CXSparse header") + SET(CXSPARSE_FOUND FALSE) +ENDIF (EXISTS ${CXSPARSE_LIB}) + +FIND_PATH(CXSPARSE_INCLUDE NAMES cs.h PATHS ${CXSPARSE_SEARCH_HEADERS}) +IF (EXISTS ${CXSPARSE_INCLUDE}) + MESSAGE("-- Found CXSparse header in: ${CXSPARSE_INCLUDE}") +ELSE (EXISTS ${CXSPARSE_INCLUDE}) + MESSAGE("-- Did not find CXSparse header") + SET(CXSPARSE_FOUND FALSE) +ENDIF (EXISTS ${CXSPARSE_INCLUDE}) + +IF (DEFINED CXSPARSE) + IF (CXSPARSE) + IF (NOT CXSPARSE_FOUND) + MESSAGE(FATAL_ERROR "-- CXSparse not found.") + ENDIF (NOT CXSPARSE_FOUND) + ELSE (CXSPARSE) + ADD_DEFINITIONS(-DCERES_NO_CXSPARSE) + ENDIF (CXSPARSE) +ELSE (DEFINED CXSPARSE) + IF (CXSPARSE_FOUND) + MESSAGE("-- Building with CXSparse support.") + SET(CXSPARSE ON) + ELSE (CXSPARSE_FOUND) + MESSAGE("-- Building without CXSparse.") + SET(CXSPARSE OFF) + ADD_DEFINITIONS(-DCERES_NO_CXSPARSE) + ENDIF (CXSPARSE_FOUND) +ENDIF (DEFINED CXSPARSE) + # Google Flags OPTION(GFLAGS "Enable Google Flags." @@ -352,6 +401,10 @@ INCLUDE_DIRECTORIES(${CHOLMOD_INCLUDE}) ENDIF(SUITESPARSE) +IF (CXSPARSE) + INCLUDE_DIRECTORIES(${CXSPARSE_INCLUDE}) +ENDIF(CXSPARSE) + IF (GFLAGS) INCLUDE_DIRECTORIES(${GFLAGS_INCLUDE}) ENDIF (GFLAGS)
diff --git a/include/ceres/solver.h b/include/ceres/solver.h index be81768..701f0b2 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h
@@ -67,11 +67,18 @@ function_tolerance = 1e-6; gradient_tolerance = 1e-10; parameter_tolerance = 1e-8; -#ifndef CERES_NO_SUITESPARSE - linear_solver_type = SPARSE_NORMAL_CHOLESKY; -#else + +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) linear_solver_type = DENSE_QR; -#endif // CERES_NO_SUITESPARSE +#else + linear_solver_type = SPARSE_NORMAL_CHOLESKY; +#endif + + sparse_linear_algebra_library = SUITE_SPARSE; +#if defined(CERES_NO_SUITESPARSE) && !defined(CERES_NO_CXSPARSE) + sparse_linear_algebra_library = CX_SPARSE; +#endif + preconditioner_type = JACOBI; num_linear_solver_threads = 1; num_eliminate_blocks = 0; @@ -154,6 +161,12 @@ // 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; + // 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/include/ceres/types.h b/include/ceres/types.h index a30c790..61a9a0d 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h
@@ -110,6 +110,15 @@ CLUSTER_TRIDIAGONAL }; +enum SparseLinearAlgebraLibraryType { + // High performance sparse Cholesky factorization and approximate + // minimum degree ordering. + SUITE_SPARSE, + + // A lightweight replacment for SuiteSparse. + CX_SPARSE +}; + enum LinearSolverTerminationType { // Termination criterion was met. For factorization based solvers // the tolerance is assumed to be zero. Any user provided values are @@ -246,6 +255,8 @@ const char* LinearSolverTypeToString(LinearSolverType type); const char* PreconditionerTypeToString(PreconditionerType type); +const char* SparseLinearAlgebraLibraryTypeToString( + SparseLinearAlgebraLibraryType type); const char* LinearSolverTerminationTypeToString( LinearSolverTerminationType type); const char* OrderingTypeToString(OrderingType type);
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 2793ee5..1b8b0d0 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -130,6 +130,10 @@ ENDIF (EXISTS ${BLAS_LIB}) ENDIF (SUITESPARSE) +IF (CXSPARSE) + LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CXSPARSE_LIB}) +ENDIF (CXSPARSE) + IF (OPENMP_FOUND) LIST(APPEND CERES_LIBRARY_DEPENDENCIES gomp) ENDIF (OPENMP_FOUND) @@ -191,13 +195,15 @@ CERES_TEST(schur_ordering) CERES_TEST(solver_impl) CERES_TEST(symmetric_linear_solver) - IF (GFLAGS) - CERES_TEST(system) - ENDIF (GFLAGS) CERES_TEST(triplet_sparse_matrix) CERES_TEST(unsymmetric_linear_solver) CERES_TEST(visibility) IF (GFLAGS) CERES_TEST(visibility_based_preconditioner) ENDIF (GFLAGS) + + # Large end to end test for the entire solver. + IF (GFLAGS) + CERES_TEST(system) + ENDIF (GFLAGS) ENDIF (BUILD_TESTING)
diff --git a/internal/ceres/linear_solver.cc b/internal/ceres/linear_solver.cc index 7221a74..c08d1f2 100644 --- a/internal/ceres/linear_solver.cc +++ b/internal/ceres/linear_solver.cc
@@ -50,22 +50,24 @@ return new CgnrSolver(options); case SPARSE_NORMAL_CHOLESKY: -#ifndef CERES_NO_SUITESPARSE - return new SparseNormalCholeskySolver(options); -#else +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) LOG(WARNING) << "SPARSE_NORMAL_CHOLESKY is not available. Please " - << "build Ceres with SuiteSparse. Returning NULL."; + << "build Ceres with SuiteSparse or CXSparse. " + << "Returning NULL."; return NULL; -#endif // CERES_NO_SUITESPARSE +#else + return new SparseNormalCholeskySolver(options); +#endif case SPARSE_SCHUR: -#ifndef CERES_NO_SUITESPARSE - return new SparseSchurComplementSolver(options); -#else +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) LOG(WARNING) << "SPARSE_SCHUR is not available. Please " - << "build Ceres with SuiteSparse. Returning NULL."; + << "build Ceres with SuiteSparse or CXSparse. " + << "Returning NULL."; return NULL; -#endif // CERES_NO_SUITESPARSE +#else + return new SparseSchurComplementSolver(options); +#endif case DENSE_SCHUR: return new DenseSchurComplementSolver(options);
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index 299a541..853abcb 100644 --- a/internal/ceres/linear_solver.h +++ b/internal/ceres/linear_solver.h
@@ -71,6 +71,7 @@ Options() : type(SPARSE_NORMAL_CHOLESKY), preconditioner_type(JACOBI), + sparse_linear_algebra_library(SUITE_SPARSE), min_num_iterations(1), max_num_iterations(1), num_threads(1), @@ -85,6 +86,8 @@ PreconditionerType preconditioner_type; + SparseLinearAlgebraLibraryType sparse_linear_algebra_library; + // Number of internal iterations that the solver uses. This // parameter only makes sense for iterative solvers like CG. int min_num_iterations;
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index 4db7c31..b30d6ed 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc
@@ -32,6 +32,11 @@ #include <ctime> #include <set> #include <vector> + +#ifndef CERES_NO_CXSPARSE +#include "cs.h" +#endif // CERES_NO_CXSPARSE + #include "Eigen/Dense" #include "ceres/block_random_access_dense_matrix.h" #include "ceres/block_random_access_matrix.h" @@ -48,6 +53,7 @@ #include "ceres/internal/scoped_ptr.h" #include "ceres/types.h" + namespace ceres { namespace internal { @@ -139,18 +145,22 @@ return true; } -#ifndef CERES_NO_SUITESPARSE + SparseSchurComplementSolver::SparseSchurComplementSolver( const LinearSolver::Options& options) - : SchurComplementSolver(options), - symbolic_factor_(NULL) { + : SchurComplementSolver(options) { +#ifndef CERES_NO_SUITESPARSE + symbolic_factor_ = NULL; +#endif // CERES_NO_SUITESPARSE } SparseSchurComplementSolver::~SparseSchurComplementSolver() { +#ifndef CERES_NO_SUITESPARSE if (symbolic_factor_ != NULL) { ss_.Free(symbolic_factor_); symbolic_factor_ = NULL; } +#endif // CERES_NO_SUITESPARSE } // Determine the non-zero blocks in the Schur Complement matrix, and @@ -224,10 +234,28 @@ set_rhs(new double[lhs()->num_rows()]); } +bool SparseSchurComplementSolver::SolveReducedLinearSystem(double* solution) { + switch (options().sparse_linear_algebra_library) { + case SUITE_SPARSE: + return SolveReducedLinearSystemUsingSuiteSparse(solution); + case CX_SPARSE: + return SolveReducedLinearSystemUsingCXSparse(solution); + default: + LOG(FATAL) << "Unknown sparse linear algebra library : " + << options().sparse_linear_algebra_library; + } + + LOG(FATAL) << "Unknown sparse linear algebra library : " + << options().sparse_linear_algebra_library; + return false; +} + +#ifndef CERES_NO_SUITESPARSE // Solve the system Sx = r, assuming that the matrix S is stored in a // BlockRandomAccessSparseMatrix. The linear system is solved using // CHOLMOD's sparse cholesky factorization routines. -bool SparseSchurComplementSolver::SolveReducedLinearSystem(double* solution) { +bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse( + double* solution) { // Extract the TripletSparseMatrix that is used for actually storing S. TripletSparseMatrix* tsm = const_cast<TripletSparseMatrix*>( @@ -272,7 +300,58 @@ ss_.Free(cholmod_solution); return true; } +#else +bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse( + double* solution) { + LOG(FATAL) << "No SuiteSparse support in Ceres."; + return false; +} #endif // CERES_NO_SUITESPARSE +#ifndef CERES_NO_CXSPARSE +// Solve the system Sx = r, assuming that the matrix S is stored in a +// BlockRandomAccessSparseMatrix. The linear system is solved using +// CXSparse's sparse cholesky factorization routines. +bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse( + double* solution) { + // Extract the TripletSparseMatrix that is used for actually storing S. + TripletSparseMatrix* tsm = + const_cast<TripletSparseMatrix*>( + down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix()); + + const int num_rows = tsm->num_rows(); + + // The case where there are no f blocks, and the system is block + // diagonal. + if (num_rows == 0) { + return true; + } + + cs_di_sparse tsm_wrapper; + tsm_wrapper.nzmax = tsm->num_nonzeros(); + tsm_wrapper.m = num_rows; + tsm_wrapper.n = num_rows; + tsm_wrapper.p = tsm->mutable_cols(); + tsm_wrapper.i = tsm->mutable_rows(); + tsm_wrapper.x = tsm->mutable_values(); + tsm_wrapper.nz = tsm->num_nonzeros(); + + cs_di_sparse* lhs = cs_compress(&tsm_wrapper); + VectorRef(solution, num_rows) = ConstVectorRef(rhs(), num_rows); + + // It maybe worth caching the ordering here, but for now we are + // going to go with the simple cholsol based implementation. + int ok = cs_di_cholsol(1, lhs, solution); + cs_free(lhs); + return ok; +} +#else +bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse( + double* solution) { + LOG(FATAL) << "No CXSparse support in Ceres."; + return false; +} +#endif // CERES_NO_CXPARSE + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h index 7cd019a..b788ce0 100644 --- a/internal/ceres/schur_complement_solver.h +++ b/internal/ceres/schur_complement_solver.h
@@ -145,7 +145,7 @@ DISALLOW_COPY_AND_ASSIGN(DenseSchurComplementSolver); }; -#ifndef CERES_NO_SUITESPARSE + // Sparse Cholesky factorization based solver. class SparseSchurComplementSolver : public SchurComplementSolver { public: @@ -155,26 +155,17 @@ private: virtual void InitStorage(const CompressedRowBlockStructure* bs); virtual bool SolveReducedLinearSystem(double* solution); + bool SolveReducedLinearSystemUsingSuiteSparse(double* solution); + bool SolveReducedLinearSystemUsingCXSparse(double* solution); - +#ifndef CERES_NO_SUITESPARSE SuiteSparse ss_; // Symbolic factorization of the reduced linear system. Precomputed // once and reused in subsequent calls. cholmod_factor* symbolic_factor_; +#endif // CERES_NO_SUITESPARSE DISALLOW_COPY_AND_ASSIGN(SparseSchurComplementSolver); }; -#else // CERES_NO_SUITESPARSE -class SparseSchurComplementSolver : public SchurComplementSolver { - public: - explicit SparseSchurComplementSolver(const LinearSolver::Options& options) - : SchurComplementSolver(options) { - LOG(FATAL) << "SPARSE_SCHUR is not available. Please " - "build Ceres with SuiteSparse."; - } - - virtual ~SparseSchurComplementSolver() {} -}; -#endif // CERES_NO_SUITESPARSE } // namespace internal } // namespace ceres
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc index 5e0cc40..ea5c346 100644 --- a/internal/ceres/schur_complement_solver_test.cc +++ b/internal/ceres/schur_complement_solver_test.cc
@@ -92,13 +92,17 @@ sol_d.get()); } - void ComputeAndCompareSolutions(int problem_id, - bool regularization, - ceres::LinearSolverType linear_solver_type) { + void ComputeAndCompareSolutions( + int problem_id, + bool regularization, + ceres::LinearSolverType linear_solver_type, + ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library) { SetUpFromProblemId(problem_id); LinearSolver::Options options; options.num_eliminate_blocks = num_eliminate_blocks; options.type = linear_solver_type; + options.sparse_linear_algebra_library = sparse_linear_algebra_library; + scoped_ptr<LinearSolver> solver(LinearSolver::Create(options)); LinearSolver::PerSolveOptions per_solve_options; @@ -133,21 +137,30 @@ }; #ifndef CERES_NO_SUITESPARSE - -TEST_F(SchurComplementSolverTest, SparseSchur) { - ComputeAndCompareSolutions(2, false, SPARSE_SCHUR); - ComputeAndCompareSolutions(3, false, SPARSE_SCHUR); - ComputeAndCompareSolutions(2, true, SPARSE_SCHUR); - ComputeAndCompareSolutions(3, true, SPARSE_SCHUR); +TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparse) { + ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE); + ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE); + ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE); + ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE); } - #endif // CERES_NO_SUITESPARSE +#ifndef CERES_NO_CXSPARSE +TEST_F(SchurComplementSolverTest, SparseSchurWithCXSparse) { + ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, CX_SPARSE); + ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, CX_SPARSE); + ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, CX_SPARSE); + ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, CX_SPARSE); +} +#endif // CERES_NO_CXSPARSE + TEST_F(SchurComplementSolverTest, DenseSchur) { - ComputeAndCompareSolutions(2, false, DENSE_SCHUR); - ComputeAndCompareSolutions(3, false, DENSE_SCHUR); - ComputeAndCompareSolutions(2, true, DENSE_SCHUR); - ComputeAndCompareSolutions(3, true, DENSE_SCHUR); + // The sparse linear algebra library type is ignored for + // DENSE_SCHUR. + ComputeAndCompareSolutions(2, false, DENSE_SCHUR, SUITE_SPARSE); + ComputeAndCompareSolutions(3, false, DENSE_SCHUR, SUITE_SPARSE); + ComputeAndCompareSolutions(2, true, DENSE_SCHUR, SUITE_SPARSE); + ComputeAndCompareSolutions(3, true, DENSE_SCHUR, SUITE_SPARSE); } } // namespace internal
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index 2ec3dee..6f8fa17 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc
@@ -432,12 +432,23 @@ LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options, string* error) { #ifdef CERES_NO_SUITESPARSE - if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY) { - *error = "Can't use SPARSE_NORMAL_CHOLESKY because SuiteSparse was not " - "enabled when Ceres was built."; + if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY && + options->sparse_linear_algebra_library == SUITE_SPARSE) { + *error = "Can't use SPARSE_NORMAL_CHOLESKY with SUITESPARSE because " + "SuiteSparse was not enabled when Ceres was built."; return NULL; } -#endif // CERES_NO_SUITESPARSE +#endif + +#ifdef CERES_NO_CXSPARSE + if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY && + options->sparse_linear_algebra_library == CX_SPARSE) { + *error = "Can't use SPARSE_NORMAL_CHOLESKY with CXSPARSE because " + "CXSparse was not enabled when Ceres was built."; + return NULL; + } +#endif + if (options->linear_solver_max_num_iterations <= 0) { *error = "Solver::Options::linear_solver_max_num_iterations is 0."; @@ -461,23 +472,25 @@ options->linear_solver_max_num_iterations; linear_solver_options.type = options->linear_solver_type; linear_solver_options.preconditioner_type = options->preconditioner_type; + linear_solver_options.sparse_linear_algebra_library = + options->sparse_linear_algebra_library; #ifdef CERES_NO_SUITESPARSE if (linear_solver_options.preconditioner_type == SCHUR_JACOBI) { *error = "SCHUR_JACOBI preconditioner not suppored. Please build Ceres " - "with SuiteSparse support"; + "with SuiteSparse support."; return NULL; } if (linear_solver_options.preconditioner_type == CLUSTER_JACOBI) { *error = "CLUSTER_JACOBI preconditioner not suppored. Please build Ceres " - "with SuiteSparse support"; + "with SuiteSparse support."; return NULL; } if (linear_solver_options.preconditioner_type == CLUSTER_TRIDIAGONAL) { *error = "CLUSTER_TRIDIAGONAL preconditioner not suppored. Please build " - "Ceres with SuiteSparse support"; + "Ceres with SuiteSparse support."; return NULL; } #endif @@ -488,23 +501,23 @@ if ((linear_solver_options.num_eliminate_blocks == 0) && IsSchurType(linear_solver_options.type)) { -#ifndef CERES_NO_SUITESPARSE +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) + LOG(INFO) << "No elimination block remaining switching to DENSE_QR."; + linear_solver_options.type = DENSE_QR; +#else LOG(INFO) << "No elimination block remaining " << "switching to SPARSE_NORMAL_CHOLESKY."; linear_solver_options.type = SPARSE_NORMAL_CHOLESKY; -#else - LOG(INFO) << "No elimination block remaining switching to DENSE_QR."; - linear_solver_options.type = DENSE_QR; -#endif // CERES_NO_SUITESPARSE +#endif } -#ifdef CERES_NO_SUITESPARSE +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) if (linear_solver_options.type == SPARSE_SCHUR) { - *error = "Can't use SPARSE_SCHUR because SuiteSparse was not " - "enabled when Ceres was built."; + *error = "Can't use SPARSE_SCHUR because neither SuiteSparse nor" + "CXSparse was enabled when Ceres was compiled."; return NULL; } -#endif // CERES_NO_SUITESPARSE +#endif // The matrix used for storing the dense Schur complement has a // single lock guarding the whole matrix. Running the
diff --git a/internal/ceres/solver_impl_test.cc b/internal/ceres/solver_impl_test.cc index 2abca63..20190f1 100644 --- a/internal/ceres/solver_impl_test.cc +++ b/internal/ceres/solver_impl_test.cc
@@ -432,14 +432,14 @@ EXPECT_EQ(parameter_blocks[2]->user_state(), &y); } -#ifdef CERES_NO_SUITESPARSE +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) TEST(SolverImpl, CreateLinearSolverNoSuiteSparse) { Solver::Options options; options.linear_solver_type = SPARSE_NORMAL_CHOLESKY; string error; EXPECT_FALSE(SolverImpl::CreateLinearSolver(&options, &error)); } -#endif // CERES_NO_SUITESPARSE +#endif TEST(SolverImpl, CreateLinearSolverNegativeMaxNumIterations) { Solver::Options options; @@ -477,11 +477,12 @@ scoped_ptr<LinearSolver> solver( SolverImpl::CreateLinearSolver(&options, &error)); EXPECT_TRUE(solver != NULL); -#ifndef CERES_NO_SUITESPARSE - EXPECT_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY); -#else + +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) EXPECT_EQ(options.linear_solver_type, DENSE_QR); -#endif // CERES_NO_SUITESPARSE +#else + EXPECT_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY); +#endif } TEST(SolverImpl, CreateLinearSolverDenseSchurMultipleThreads) { @@ -508,10 +509,19 @@ #ifndef CERES_NO_SUITESPARSE options.linear_solver_type = SPARSE_NORMAL_CHOLESKY; + options.sparse_linear_algebra_library = SUITE_SPARSE; solver.reset(SolverImpl::CreateLinearSolver(&options, &error)); EXPECT_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY); EXPECT_TRUE(solver.get() != NULL); -#endif // CERES_NO_SUITESPARSE +#endif + +#ifndef CERES_NO_CXSPARSE + options.linear_solver_type = SPARSE_NORMAL_CHOLESKY; + options.sparse_linear_algebra_library = CX_SPARSE; + solver.reset(SolverImpl::CreateLinearSolver(&options, &error)); + EXPECT_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY); + EXPECT_TRUE(solver.get() != NULL); +#endif options.linear_solver_type = DENSE_SCHUR; options.num_eliminate_blocks = 2; @@ -521,13 +531,14 @@ options.linear_solver_type = SPARSE_SCHUR; options.num_eliminate_blocks = 2; -#ifndef CERES_NO_SUITESPARSE solver.reset(SolverImpl::CreateLinearSolver(&options, &error)); + +#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) + EXPECT_TRUE(SolverImpl::CreateLinearSolver(&options, &error) == NULL); +#else EXPECT_TRUE(solver.get() != NULL); EXPECT_EQ(options.linear_solver_type, SPARSE_SCHUR); -#else // CERES_NO_SUITESPARSE - EXPECT_TRUE(SolverImpl::CreateLinearSolver(&options, &error) == NULL); -#endif // CERES_NO_SUITESPARSE +#endif options.linear_solver_type = ITERATIVE_SCHUR; options.num_eliminate_blocks = 2;
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc index d4c2f3d..5a082de 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -28,13 +28,16 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) -#ifndef CERES_NO_SUITESPARSE - #include "ceres/sparse_normal_cholesky_solver.h" #include <algorithm> #include <cstring> #include <ctime> + +#ifndef CERES_NO_CXSPARSE +#include "cs.h" +#endif + #include "ceres/compressed_row_sparse_matrix.h" #include "ceres/linear_solver.h" #include "ceres/suitesparse.h" @@ -43,18 +46,26 @@ #include "ceres/internal/scoped_ptr.h" #include "ceres/types.h" + + namespace ceres { namespace internal { SparseNormalCholeskySolver::SparseNormalCholeskySolver( const LinearSolver::Options& options) - : options_(options), symbolic_factor_(NULL) {} + : options_(options) { +#ifndef CERES_NO_SUITESPARSE + symbolic_factor_ = NULL; +#endif +} SparseNormalCholeskySolver::~SparseNormalCholeskySolver() { +#ifndef CERES_NO_SUITESPARSE if (symbolic_factor_ != NULL) { ss_.Free(symbolic_factor_); symbolic_factor_ = NULL; } +#endif } LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl( @@ -62,6 +73,94 @@ const double* b, const LinearSolver::PerSolveOptions& per_solve_options, double * x) { + switch (options_.sparse_linear_algebra_library) { + case SUITE_SPARSE: + return SolveImplUsingSuiteSparse(A, b, per_solve_options, x); + case CX_SPARSE: + return SolveImplUsingCXSparse(A, b, per_solve_options, x); + default: + LOG(FATAL) << "Unknown sparse linear algebra library : " + << options_.sparse_linear_algebra_library; + } + + LOG(FATAL) << "Unknown sparse linear algebra library : " + << options_.sparse_linear_algebra_library; + return LinearSolver::Summary(); +} + +#ifndef CERES_NO_CXSPARSE +LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double * x) { + LinearSolver::Summary summary; + summary.num_iterations = 1; + const int num_cols = A->num_cols(); + Vector Atb = Vector::Zero(num_cols); + A->LeftMultiply(b, Atb.data()); + + 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. + CompressedRowSparseMatrix D(per_solve_options.D, num_cols); + A->AppendRows(D); + } + + VectorRef(x, num_cols).setZero(); + + // Wrap the augmented Jacobian in a compressed sparse column matrix. + cs_di At; + At.m = A->num_cols(); + At.n = A->num_rows(); + At.nz = -1; + At.nzmax = A->num_nonzeros(); + At.p = A->mutable_rows(); + At.i = A->mutable_cols(); + At.x = A->mutable_values(); + + // 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. + cs_di* A2 = cs_transpose(&At, 1); + cs_di* AtA = cs_multiply(&At,A2); + + cs_free(A2); + if (per_solve_options.D != NULL) { + A->DeleteRows(num_cols); + } + + // This recomputes the symbolic factorization every time it is + // invoked. It will perhaps be worth it to cache the symbolic + // factorization the way we do for SuiteSparse. + if (cs_cholsol(1, AtA, Atb.data())) { + VectorRef(x, Atb.rows()) = Atb; + summary.termination_type = TOLERANCE; + } + + cs_free(AtA); + return summary; +} +#else +LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double * x) { + LOG(FATAL) << "No CXSparse support in Ceres."; +} +#endif + +#ifndef CERES_NO_SUITESPARSE +LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double * x) { const time_t start_time = time(NULL); const int num_cols = A->num_cols(); @@ -117,8 +216,15 @@ << " cleanup: " << cleanup_time - solve_time; return summary; } +#else +LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double * x) { + LOG(FATAL) << "No SuiteSparse support in Ceres."; +} +#endif } // namespace internal } // namespace ceres - -#endif // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h index ce1d6d2..5445559 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.h +++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -35,11 +35,11 @@ #define CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_ #ifndef CERES_NO_SUITESPARSE - #include "cholmod.h" -#include "cholmod_core.h" -#include "ceres/linear_solver.h" #include "ceres/suitesparse.h" +#endif // CERES_NO_SUITESPARSE + +#include "ceres/linear_solver.h" #include "ceres/internal/macros.h" namespace ceres { @@ -61,17 +61,31 @@ const LinearSolver::PerSolveOptions& options, double* x); - const LinearSolver::Options options_; - SuiteSparse ss_; + LinearSolver::Summary SolveImplUsingSuiteSparse( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& options, + double* x); + // Crashes if CSparse is not installed. + LinearSolver::Summary SolveImplUsingCXSparse( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& options, + double* x); + +#ifndef CERES_NO_SUITESPARSE + SuiteSparse ss_; // Cached factorization cholmod_factor* symbolic_factor_; +#endif // CERES_NO_SUITESPARSE + + + const LinearSolver::Options options_; DISALLOW_COPY_AND_ASSIGN(SparseNormalCholeskySolver); }; } // namespace internal } // namespace ceres -#endif // CERES_NO_SUITESPARSE - #endif // CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
diff --git a/internal/ceres/system_test.cc b/internal/ceres/system_test.cc index a9a6747..8e2f5eb 100644 --- a/internal/ceres/system_test.cc +++ b/internal/ceres/system_test.cc
@@ -62,32 +62,39 @@ // Struct used for configuring the solver. struct SolverConfig { SolverConfig(LinearSolverType linear_solver_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library, OrderingType ordering_type) : linear_solver_type(linear_solver_type), + sparse_linear_algebra_library(sparse_linear_algebra_library), ordering_type(ordering_type), preconditioner_type(IDENTITY), num_threads(1) { } SolverConfig(LinearSolverType linear_solver_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library, OrderingType ordering_type, PreconditionerType preconditioner_type, int num_threads) : linear_solver_type(linear_solver_type), + sparse_linear_algebra_library(sparse_linear_algebra_library), ordering_type(ordering_type), preconditioner_type(preconditioner_type), num_threads(num_threads) { } string ToString() const { - return StringPrintf("(%s, %s, %s, %d)", - LinearSolverTypeToString(linear_solver_type), - OrderingTypeToString(ordering_type), - PreconditionerTypeToString(preconditioner_type), - num_threads); + return StringPrintf( + "(%s, %s, %s, %s, %d)", + LinearSolverTypeToString(linear_solver_type), + SparseLinearAlgebraLibraryTypeToString(sparse_linear_algebra_library), + OrderingTypeToString(ordering_type), + PreconditionerTypeToString(preconditioner_type), + num_threads); } LinearSolverType linear_solver_type; + SparseLinearAlgebraLibraryType sparse_linear_algebra_library; OrderingType ordering_type; PreconditionerType preconditioner_type; int num_threads; @@ -122,6 +129,8 @@ Solver::Options& options = *(system_test_problem->mutable_solver_options()); options.linear_solver_type = config.linear_solver_type; + options.sparse_linear_algebra_library = + config.sparse_linear_algebra_library; options.ordering_type = config.ordering_type; options.preconditioner_type = config.preconditioner_type; options.num_threads = config.num_threads; @@ -266,19 +275,31 @@ }; TEST(SystemTest, PowellsFunction) { - vector<SolverConfig> configurations; - configurations.push_back(SolverConfig(DENSE_QR, NATURAL)); - configurations.push_back(SolverConfig(DENSE_SCHUR, SCHUR)); + vector<SolverConfig> configs; +#define CONFIGURE(linear_solver, sparse_linear_algebra_library, ordering) \ + configs.push_back(SolverConfig(linear_solver, \ + sparse_linear_algebra_library, \ + ordering)) + + CONFIGURE(DENSE_QR, SUITE_SPARSE, NATURAL); + CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, SCHUR); #ifndef CERES_NO_SUITESPARSE - configurations.push_back(SolverConfig(SPARSE_NORMAL_CHOLESKY, NATURAL)); - configurations.push_back(SolverConfig(SPARSE_SCHUR, SCHUR)); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, NATURAL); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, SCHUR); #endif // CERES_NO_SUITESPARSE - configurations.push_back(SolverConfig(ITERATIVE_SCHUR, SCHUR)); +#ifndef CERES_NO_CXSPARSE + CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, NATURAL); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, SCHUR); +#endif // CERES_NO_CXSPARSE + + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, SCHUR); + +#undef CONFIGURE + const double kMaxAbsoluteDifference = 1e-8; - RunSolversAndCheckTheyMatch<PowellsFunction>(configurations, - kMaxAbsoluteDifference); + RunSolversAndCheckTheyMatch<PowellsFunction>(configs, kMaxAbsoluteDifference); } // This class implements the SystemTestProblem interface and provides @@ -462,36 +483,50 @@ TEST(SystemTest, BundleAdjustmentProblem) { vector<SolverConfig> configs; -#define CONFIGURE(a, b, c, d) configs.push_back(SolverConfig(a, b, c, d)) +#define CONFIGURE(linear_solver, sparse_linear_algebra_library, ordering, preconditioner, threads) \ + configs.push_back(SolverConfig(linear_solver, \ + sparse_linear_algebra_library, \ + ordering, \ + preconditioner, \ + threads)) #ifndef CERES_NO_SUITESPARSE - CONFIGURE(SPARSE_NORMAL_CHOLESKY, NATURAL, IDENTITY, 1); - CONFIGURE(SPARSE_NORMAL_CHOLESKY, USER, IDENTITY, 1); - CONFIGURE(SPARSE_NORMAL_CHOLESKY, SCHUR, IDENTITY, 1); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, NATURAL, IDENTITY, 1); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, USER, IDENTITY, 1); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, SCHUR, IDENTITY, 1); - CONFIGURE(SPARSE_SCHUR, USER, IDENTITY, 1); - CONFIGURE(SPARSE_SCHUR, SCHUR, IDENTITY, 1); + CONFIGURE(SPARSE_SCHUR, SUITE_SPARSE, USER, IDENTITY, 1); + CONFIGURE(SPARSE_SCHUR, SUITE_SPARSE, SCHUR, IDENTITY, 1); #endif // CERES_NO_SUITESPARSE - CONFIGURE(DENSE_SCHUR, USER, IDENTITY, 1); - CONFIGURE(DENSE_SCHUR, SCHUR, IDENTITY, 1); +#ifndef CERES_NO_CXSPARSE + CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, NATURAL, IDENTITY, 1); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, USER, IDENTITY, 1); + CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, SCHUR, IDENTITY, 1); - CONFIGURE(CGNR, USER, JACOBI, 1); + CONFIGURE(SPARSE_SCHUR, CX_SPARSE, USER, IDENTITY, 1); + CONFIGURE(SPARSE_SCHUR, CX_SPARSE, SCHUR, IDENTITY, 1); +#endif // CERES_NO_CXSPARSE - CONFIGURE(ITERATIVE_SCHUR, USER, JACOBI, 1); + CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, USER, IDENTITY, 1); + CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, SCHUR, IDENTITY, 1); + + CONFIGURE(CGNR, SUITE_SPARSE, USER, JACOBI, 1); + + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, USER, JACOBI, 1); #ifndef CERES_NO_SUITESPARSE - CONFIGURE(ITERATIVE_SCHUR, USER, SCHUR_JACOBI, 1); - CONFIGURE(ITERATIVE_SCHUR, USER, CLUSTER_JACOBI, 1); - CONFIGURE(ITERATIVE_SCHUR, USER, CLUSTER_TRIDIAGONAL, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, USER, SCHUR_JACOBI, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, USER, CLUSTER_JACOBI, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, USER, CLUSTER_TRIDIAGONAL, 1); #endif // CERES_NO_SUITESPARSE - CONFIGURE(ITERATIVE_SCHUR, SCHUR, JACOBI, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, SCHUR, JACOBI, 1); #ifndef CERES_NO_SUITESPARSE - CONFIGURE(ITERATIVE_SCHUR, SCHUR, SCHUR_JACOBI, 1); - CONFIGURE(ITERATIVE_SCHUR, SCHUR, CLUSTER_JACOBI, 1); - CONFIGURE(ITERATIVE_SCHUR, SCHUR, CLUSTER_TRIDIAGONAL, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, SCHUR, SCHUR_JACOBI, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, SCHUR, CLUSTER_JACOBI, 1); + CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, SCHUR, CLUSTER_TRIDIAGONAL, 1); #endif // CERES_NO_SUITESPARSE #undef CONFIGURE
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 860f8a4..c6c4dea 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc
@@ -61,6 +61,16 @@ } } +const char* SparseLinearAlgebraLibraryTypeToString( + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type) { + switch (sparse_linear_algebra_library_type) { + CASESTR(SUITE_SPARSE); + CASESTR(CX_SPARSE); + default: + return "UNKNOWN"; + } +} + const char* OrderingTypeToString(OrderingType ordering_type) { switch (ordering_type) { CASESTR(NATURAL);
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc index be91056..a7ee23e 100644 --- a/internal/ceres/unsymmetric_linear_solver_test.cc +++ b/internal/ceres/unsymmetric_linear_solver_test.cc
@@ -52,109 +52,82 @@ A_.reset(down_cast<TripletSparseMatrix*>(problem->A.release())); b_.reset(problem->b.release()); D_.reset(problem->D.release()); - sol1_.reset(problem->x.release()); - sol2_.reset(problem->x_D.release()); - x_.reset(new double[A_->num_cols()]); + sol_unregularized_.reset(problem->x.release()); + sol_regularized_.reset(problem->x_D.release()); } - void TestSolver(LinearSolverType linear_solver_type) { + void TestSolver( + LinearSolverType linear_solver_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library) { LinearSolver::Options options; options.type = linear_solver_type; + options.sparse_linear_algebra_library = sparse_linear_algebra_library; scoped_ptr<LinearSolver> solver(LinearSolver::Create(options)); LinearSolver::PerSolveOptions per_solve_options; + LinearSolver::Summary unregularized_solve_summary; + LinearSolver::Summary regularized_solve_summary; + Vector x_unregularized(A_->num_cols()); + Vector x_regularized(A_->num_cols()); - // Unregularized - LinearSolver::Summary summary = - solver->Solve(A_.get(), b_.get(), per_solve_options, x_.get()); + scoped_ptr<SparseMatrix> transformed_A; - EXPECT_EQ(summary.termination_type, TOLERANCE); - - for (int i = 0; i < A_->num_cols(); ++i) { - EXPECT_NEAR(sol1_[i], x_[i], 1e-8); + if (linear_solver_type == DENSE_QR) { + transformed_A.reset(new DenseSparseMatrix(*A_)); + } else if (linear_solver_type == SPARSE_NORMAL_CHOLESKY) { + transformed_A.reset(new CompressedRowSparseMatrix(*A_)); + } else { + LOG(FATAL) << "Unknown linear solver : " << linear_solver_type; } + // Unregularized + unregularized_solve_summary = + solver->Solve(transformed_A.get(), + b_.get(), + per_solve_options, + x_unregularized.data()); // Regularized solution per_solve_options.D = D_.get(); - summary = solver->Solve(A_.get(), b_.get(), per_solve_options, x_.get()); + regularized_solve_summary = + solver->Solve(transformed_A.get(), + b_.get(), + per_solve_options, + x_regularized.data()); - EXPECT_EQ(summary.termination_type, TOLERANCE); + EXPECT_EQ(unregularized_solve_summary.termination_type, TOLERANCE); for (int i = 0; i < A_->num_cols(); ++i) { - EXPECT_NEAR(sol2_[i], x_[i], 1e-8); + EXPECT_NEAR(sol_unregularized_[i], x_unregularized[i], 1e-8); + } + + EXPECT_EQ(regularized_solve_summary.termination_type, TOLERANCE); + for (int i = 0; i < A_->num_cols(); ++i) { + EXPECT_NEAR(sol_regularized_[i], x_regularized[i], 1e-8); } } scoped_ptr<TripletSparseMatrix> A_; scoped_array<double> b_; scoped_array<double> D_; - scoped_array<double> sol1_; - scoped_array<double> sol2_; - - scoped_array<double> x_; + scoped_array<double> sol_unregularized_; + scoped_array<double> sol_regularized_; }; -// TODO(keir): Reduce duplication. TEST_F(UnsymmetricLinearSolverTest, DenseQR) { - LinearSolver::Options options; - options.type = DENSE_QR; - scoped_ptr<LinearSolver> solver(LinearSolver::Create(options)); - - LinearSolver::PerSolveOptions per_solve_options; - DenseSparseMatrix A(*A_); - - // Unregularized - LinearSolver::Summary summary = - solver->Solve(&A, b_.get(), per_solve_options, x_.get()); - - EXPECT_EQ(summary.termination_type, TOLERANCE); - for (int i = 0; i < A_->num_cols(); ++i) { - EXPECT_NEAR(sol1_[i], x_[i], 1e-8); - } - - VectorRef x(x_.get(), A_->num_cols()); - VectorRef b(b_.get(), A_->num_rows()); - Vector r = A.matrix()*x - b; - LOG(INFO) << "r = A*x - b: \n" << r; - - // Regularized solution - per_solve_options.D = D_.get(); - summary = solver->Solve(&A, b_.get(), per_solve_options, x_.get()); - - EXPECT_EQ(summary.termination_type, TOLERANCE); - for (int i = 0; i < A_->num_cols(); ++i) { - EXPECT_NEAR(sol2_[i], x_[i], 1e-8); - } + TestSolver(DENSE_QR, SUITE_SPARSE); } #ifndef CERES_NO_SUITESPARSE -TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholesky) { - LinearSolver::Options options; - options.type = SPARSE_NORMAL_CHOLESKY; - scoped_ptr<LinearSolver>solver(LinearSolver::Create(options)); - - LinearSolver::PerSolveOptions per_solve_options; - CompressedRowSparseMatrix A(*A_); - - // Unregularized - LinearSolver::Summary summary = - solver->Solve(&A, b_.get(), per_solve_options, x_.get()); - - EXPECT_EQ(summary.termination_type, TOLERANCE); - for (int i = 0; i < A_->num_cols(); ++i) { - EXPECT_NEAR(sol1_[i], x_[i], 1e-8); - } - - // Regularized solution - per_solve_options.D = D_.get(); - summary = solver->Solve(&A, b_.get(), per_solve_options, x_.get()); - - EXPECT_EQ(summary.termination_type, TOLERANCE); - for (int i = 0; i < A_->num_cols(); ++i) { - EXPECT_NEAR(sol2_[i], x_[i], 1e-8); - } +TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingSuiteSparse) { + TestSolver(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE); } -#endif // CERES_NO_SUITESPARSE +#endif + +#ifndef CERES_NO_CXSPARSE +TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingCXSparse) { + TestSolver(SPARSE_NORMAL_CHOLESKY, CX_SPARSE); +} +#endif } // namespace internal } // namespace ceres
diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h index 5f14842..c265a0e 100644 --- a/internal/ceres/visibility_based_preconditioner.h +++ b/internal/ceres/visibility_based_preconditioner.h
@@ -261,7 +261,7 @@ virtual void LeftMultiply(const double* x, double* y) const {} virtual int num_rows() const { return -1; } virtual int num_cols() const { return -1; } - bool Compute(const BlockSparseMatrixBase& A, const double* D) { + bool Update(const BlockSparseMatrixBase& A, const double* D) { return false; } };