Caching the symbolic Cholesky factorization when using CXSparse
Average factorization times for bundle adjustment test problem:
SuiteSparse: 0.2794 s.
CXSparse: 0.4039 s.
CXSparse cached: 0.2399 s.
CXSparse will still be slower, though, because it has to compute
the transpose and J^T * J.
Change-Id: If9cdaa3dd520bee84b56e5fd4953b56a93db6bde
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index a7c43ef..9e00b44 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -55,6 +55,10 @@
#ifndef CERES_NO_SUITESPARSE
factor_ = NULL;
#endif
+
+#ifndef CERES_NO_CXSPARSE
+ cxsparse_factor_ = NULL;
+#endif // CERES_NO_CXSPARSE
}
SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {
@@ -64,6 +68,13 @@
factor_ = NULL;
}
#endif
+
+#ifndef CERES_NO_CXSPARSE
+ if (cxsparse_factor_ != NULL) {
+ cxsparse_.Free(cxsparse_factor_);
+ cxsparse_factor_ = NULL;
+ }
+#endif // CERES_NO_CXSPARSE
}
LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
@@ -108,14 +119,7 @@
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();
+ cs_di At = cxsparse_.CreateSparseMatrixTransposeView(A);
// Compute the normal equations. J'J delta = J'f and solve them
// using a sparse Cholesky factorization. Notice that when compared
@@ -127,20 +131,23 @@
cs_di* A2 = cs_transpose(&At, 1);
cs_di* AtA = cs_multiply(&At,A2);
- cs_free(A2);
+ cxsparse_.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())) {
+ // Compute symbolic factorization if not available.
+ if (cxsparse_factor_ == NULL) {
+ cxsparse_factor_ = CHECK_NOTNULL(cxsparse_.AnalyzeCholesky(AtA));
+ }
+
+ // Solve the linear system.
+ if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) {
VectorRef(x, Atb.rows()) = Atb;
summary.termination_type = TOLERANCE;
}
- cs_free(AtA);
+ cxsparse_.Free(AtA);
return summary;
}
#else
