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/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;
}
};
