Refactor suitesparse.h/cc
1. Generalize SuiteSparse::AnalyzeCholesky and
SuiteSparse::BlockAnalyzeCholesky from just doing AMD to taking
OrderingType as an argument and using that to determine whether
AMD & Nested Dissection algorithms are used for computing the
fill-reducing ordering or a natural ordering when computing
the symbolic factorization.
2. Remove AnalyzeCholeskyWithNaturalOrdering.
3. Replace and generalize SuiteSparse::BlockAMDOrdering with
SuiteSparse::BlockOrdering which also takes OrderingType as an
argument. Same for SuiteSparse::ApproximateMinimumDegreeOrdering
and SuiteSparse::NestedDissectionOrdering by
SuiteSparse::Ordering.
4. Remove LinearSolver::Options::use_postordering and replace it
with LinearSolver::Options::ordering_type.
5. Replace Preconditioner::Options::use_postordering and replace it
with Preconditioner::Options::ordering_type.
6. Add NESDIS to OrderingType. With the above changes, the linear
solvers can now use Nested Dissection once this information
is piped through the nonlinear solver.
Change-Id: Ib8e93fbf34ae2981bf2ac54dcda9e25c7c213790
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc
index c7cf4f8..12e2ef9 100644
--- a/internal/ceres/cgnr_solver.cc
+++ b/internal/ceres/cgnr_solver.cc
@@ -80,7 +80,7 @@
options_.subset_preconditioner_start_row_block;
preconditioner_options.sparse_linear_algebra_library_type =
options_.sparse_linear_algebra_library_type;
- preconditioner_options.use_postordering = options_.use_postordering;
+ preconditioner_options.ordering_type = options_.ordering_type;
preconditioner_options.num_threads = options_.num_threads;
preconditioner_options.context = options_.context;
preconditioner_ =
diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
index f31004b..3abfb29 100644
--- a/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver.cc
@@ -254,7 +254,8 @@
const int num_cols = A->num_cols();
cholmod_sparse lhs = ss.CreateSparseMatrixTransposeView(A);
event_logger.AddEvent("Setup");
- cholmod_factor* factor = ss.AnalyzeCholesky(&lhs, &summary.message);
+ cholmod_factor* factor =
+ ss.AnalyzeCholesky(&lhs, options_.ordering_type, &summary.message);
event_logger.AddEvent("Analysis");
if (factor == nullptr) {
diff --git a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
index f3bf650..2621690 100644
--- a/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
+++ b/internal/ceres/dynamic_sparse_normal_cholesky_solver_test.cc
@@ -93,12 +93,14 @@
}
void TestSolver(
- const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type) {
+ const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+ const OrderingType ordering_type) {
LinearSolver::Options options;
options.type = SPARSE_NORMAL_CHOLESKY;
options.dynamic_sparsity = true;
options.sparse_linear_algebra_library_type =
sparse_linear_algebra_library_type;
+ options.ordering_type = ordering_type;
ContextImpl context;
options.context = &context;
TestSolver(options, nullptr);
@@ -111,20 +113,25 @@
};
#ifndef CERES_NO_SUITESPARSE
-TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparse) {
- TestSolver(SUITE_SPARSE);
+TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparseAMD) {
+ TestSolver(SUITE_SPARSE, OrderingType::AMD);
}
+#ifndef CERES_NO_METIS
+TEST_F(DynamicSparseNormalCholeskySolverTest, SuiteSparseNESDIS) {
+ TestSolver(SUITE_SPARSE, OrderingType::NESDIS);
+}
+#endif
#endif
#ifndef CERES_NO_CXSPARSE
TEST_F(DynamicSparseNormalCholeskySolverTest, CXSparse) {
- TestSolver(CX_SPARSE);
+ TestSolver(CX_SPARSE, OrderingType::AMD);
}
#endif
#ifdef CERES_USE_EIGEN_SPARSE
TEST_F(DynamicSparseNormalCholeskySolverTest, Eigen) {
- TestSolver(EIGEN_SPARSE);
+ TestSolver(EIGEN_SPARSE, OrderingType::AMD);
}
#endif // CERES_USE_EIGEN_SPARSE
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index 57c1a79..4916e41 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -72,6 +72,7 @@
// the wrong arguments.
FATAL_ERROR
};
+
inline std::ostream& operator<<(std::ostream& s,
LinearSolverTerminationType type) {
switch (type) {
@@ -105,8 +106,27 @@
AMD, // Use the Approximate Minimum Degree algorithm to re-order
// the matrix.
+
+ NESDIS, // Use the Nested Dissection algorithm to re-order the matrix.
};
+inline std::ostream& operator<<(std::ostream& s, OrderingType type) {
+ switch (type) {
+ case OrderingType::NATURAL:
+ s << "NATURAL";
+ break;
+ case OrderingType::AMD:
+ s << "AMD";
+ break;
+ case OrderingType::NESDIS:
+ s << "NESDIS";
+ break;
+ default:
+ s << "UNKNOWN OrderingType";
+ }
+ return s;
+}
+
class LinearOperator;
// Abstract base class for objects that implement algorithms for
@@ -134,9 +154,9 @@
DenseLinearAlgebraLibraryType dense_linear_algebra_library_type = EIGEN;
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type =
SUITE_SPARSE;
+ OrderingType ordering_type = OrderingType::NATURAL;
// See solver.h for information about these flags.
- bool use_postordering = false;
bool dynamic_sparsity = false;
bool use_explicit_schur_complement = false;
diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h
index fdc4d86..68b575f 100644
--- a/internal/ceres/preconditioner.h
+++ b/internal/ceres/preconditioner.h
@@ -39,6 +39,7 @@
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/export.h"
#include "ceres/linear_operator.h"
+#include "ceres/linear_solver.h"
#include "ceres/sparse_matrix.h"
#include "ceres/types.h"
@@ -54,6 +55,7 @@
VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS;
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type =
SUITE_SPARSE;
+ OrderingType ordering_type = OrderingType::NATURAL;
// When using the subset preconditioner, all row blocks starting
// from this row block are used to construct the preconditioner.
@@ -67,9 +69,6 @@
// and the preconditioner is the inverse of the matrix Q'Q.
int subset_preconditioner_start_row_block = -1;
- // See solver.h for information about these flags.
- bool use_postordering = false;
-
// If possible, how many threads the preconditioner can use.
int num_threads = 1;
diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc
index bda4545..c6765fe 100644
--- a/internal/ceres/reorder_program.cc
+++ b/internal/ceres/reorder_program.cc
@@ -123,7 +123,7 @@
// If the user did not supply a useful ordering, then just use
// regular AMD.
if (parameter_block_ordering.NumGroups() <= 1) {
- ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose, &ordering[0]);
+ ss.Ordering(block_jacobian_transpose, OrderingType::AMD, &ordering[0]);
} else {
vector<int> constraints;
for (auto* parameter_block : parameter_blocks) {
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc
index 697d54a..362a1f6 100644
--- a/internal/ceres/schur_complement_solver_test.cc
+++ b/internal/ceres/schur_complement_solver_test.cc
@@ -93,7 +93,7 @@
ceres::LinearSolverType linear_solver_type,
ceres::DenseLinearAlgebraLibraryType dense_linear_algebra_library_type,
ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
- bool use_postordering) {
+ ceres::internal::OrderingType ordering_type) {
SetUpFromProblemId(problem_id);
LinearSolver::Options options;
options.elimination_groups.push_back(num_eliminate_blocks);
@@ -104,7 +104,7 @@
dense_linear_algebra_library_type;
options.sparse_linear_algebra_library_type =
sparse_linear_algebra_library_type;
- options.use_postordering = use_postordering;
+ options.ordering_type = ordering_type;
ContextImpl context;
options.context = &context;
DetectStructure(*A->block_structure(),
@@ -150,96 +150,154 @@
// TODO(sameeragarwal): Refactor these using value parameterized tests.
// TODO(sameeragarwal): More extensive tests using random matrices.
TEST_F(SchurComplementSolverTest, DenseSchurWithEigenSmallProblem) {
- ComputeAndCompareSolutions(2, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
- ComputeAndCompareSolutions(2, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 2, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 2, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
}
TEST_F(SchurComplementSolverTest, DenseSchurWithEigenLargeProblem) {
- ComputeAndCompareSolutions(3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
- ComputeAndCompareSolutions(3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
}
TEST_F(SchurComplementSolverTest, DenseSchurWithEigenVaryingFBlockSize) {
- ComputeAndCompareSolutions(4, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 4, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
}
#ifndef CERES_NO_LAPACK
TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKSmallProblem) {
- ComputeAndCompareSolutions(2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
- ComputeAndCompareSolutions(2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL);
}
TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKLargeProblem) {
- ComputeAndCompareSolutions(3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
- ComputeAndCompareSolutions(3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, OrderingType::NATURAL);
}
#endif
#ifndef CERES_NO_SUITESPARSE
TEST_F(SchurComplementSolverTest,
- SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) {
+ SparseSchurWithSuiteSparseSmallProblemNATURAL) {
ComputeAndCompareSolutions(
- 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
- ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
-}
-
-TEST_F(SchurComplementSolverTest,
- SparseSchurWithSuiteSparseSmallProblemPostOrdering) {
- ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
- ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
-}
-
-TEST_F(SchurComplementSolverTest,
- SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) {
+ 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
ComputeAndCompareSolutions(
- 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
- ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
+ 2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
}
TEST_F(SchurComplementSolverTest,
- SparseSchurWithSuiteSparseLargeProblemPostOrdering) {
- ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
- ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
+ SparseSchurWithSuiteSparseLargeProblemNATURAL) {
+ ComputeAndCompareSolutions(
+ 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NATURAL);
}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseSmallProblemAMD) {
+ ComputeAndCompareSolutions(
+ 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD);
+ ComputeAndCompareSolutions(
+ 2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD);
+}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseLargeProblemAMD) {
+ ComputeAndCompareSolutions(
+ 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD);
+ ComputeAndCompareSolutions(
+ 3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::AMD);
+}
+
+#ifndef CERES_NO_METIS
+TEST_F(SchurComplementSolverTest,
+ SparseSchurWithSuiteSparseSmallProblemNESDIS) {
+ ComputeAndCompareSolutions(
+ 2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS);
+ ComputeAndCompareSolutions(
+ 2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS);
+}
+TEST_F(SchurComplementSolverTest,
+ SparseSchurWithSuiteSparseLargeProblemNESDIS) {
+ ComputeAndCompareSolutions(
+ 3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS);
+ ComputeAndCompareSolutions(
+ 3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, OrderingType::NESDIS);
+}
+#endif // CERES_NO_METIS
#endif // CERES_NO_SUITESPARSE
#ifndef CERES_NO_CXSPARSE
+// TODO(sameeragarwal): Extend these tests for NATURAL & NESDIS, once the linear
+// solver supports it.
TEST_F(SchurComplementSolverTest, SparseSchurWithCXSparseSmallProblem) {
- ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
- ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 2, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD);
+ ComputeAndCompareSolutions(
+ 2, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD);
}
TEST_F(SchurComplementSolverTest, SparseSchurWithCXSparseLargeProblem) {
- ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
- ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
+ ComputeAndCompareSolutions(
+ 3, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD);
+ ComputeAndCompareSolutions(
+ 3, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, OrderingType::AMD);
}
#endif // CERES_NO_CXSPARSE
#ifndef CERES_NO_ACCELERATE_SPARSE
+// TODO(sameeragarwal): Extend these tests for NATURAL & NESDIS, once the linear
+// solver supports it.
TEST_F(SchurComplementSolverTest, SparseSchurWithAccelerateSparseSmallProblem) {
ComputeAndCompareSolutions(
- 2, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true);
+ 2, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
ComputeAndCompareSolutions(
- 2, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true);
+ 2, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
}
TEST_F(SchurComplementSolverTest, SparseSchurWithAccelerateSparseLargeProblem) {
ComputeAndCompareSolutions(
- 3, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true);
+ 3, false, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
ComputeAndCompareSolutions(
- 3, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, true);
+ 3, true, SPARSE_SCHUR, EIGEN, ACCELERATE_SPARSE, OrderingType::AMD);
}
#endif // CERES_NO_ACCELERATE_SPARSE
#ifdef CERES_USE_EIGEN_SPARSE
-TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseSmallProblem) {
- ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
- ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
+TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseSmallProblemAMD) {
+ ComputeAndCompareSolutions(
+ 2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD);
+ ComputeAndCompareSolutions(
+ 2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD);
}
-TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseLargeProblem) {
- ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
- ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
+TEST_F(SchurComplementSolverTest,
+ SparseSchurWithEigenSparseSmallProblemNATURAL) {
+ ComputeAndCompareSolutions(
+ 2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL);
+}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithEigenSparseLargeProblemAMD) {
+ ComputeAndCompareSolutions(
+ 3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD);
+ ComputeAndCompareSolutions(
+ 3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::AMD);
+}
+
+TEST_F(SchurComplementSolverTest,
+ SparseSchurWithEigenSparseLargeProblemNATURAL) {
+ ComputeAndCompareSolutions(
+ 3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL);
+ ComputeAndCompareSolutions(
+ 3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, OrderingType::NATURAL);
}
#endif // CERES_USE_EIGEN_SPARSE
diff --git a/internal/ceres/sparse_cholesky.cc b/internal/ceres/sparse_cholesky.cc
index ac7c9b2..c68b833 100644
--- a/internal/ceres/sparse_cholesky.cc
+++ b/internal/ceres/sparse_cholesky.cc
@@ -44,17 +44,16 @@
std::unique_ptr<SparseCholesky> SparseCholesky::Create(
const LinearSolver::Options& options) {
- const OrderingType ordering_type =
- options.use_postordering ? OrderingType::AMD : OrderingType::NATURAL;
std::unique_ptr<SparseCholesky> sparse_cholesky;
switch (options.sparse_linear_algebra_library_type) {
case SUITE_SPARSE:
#ifndef CERES_NO_SUITESPARSE
if (options.use_mixed_precision_solves) {
- sparse_cholesky = FloatSuiteSparseCholesky::Create(ordering_type);
+ sparse_cholesky =
+ FloatSuiteSparseCholesky::Create(options.ordering_type);
} else {
- sparse_cholesky = SuiteSparseCholesky::Create(ordering_type);
+ sparse_cholesky = SuiteSparseCholesky::Create(options.ordering_type);
}
break;
#else
@@ -64,9 +63,10 @@
case EIGEN_SPARSE:
#ifdef CERES_USE_EIGEN_SPARSE
if (options.use_mixed_precision_solves) {
- sparse_cholesky = FloatEigenSparseCholesky::Create(ordering_type);
+ sparse_cholesky =
+ FloatEigenSparseCholesky::Create(options.ordering_type);
} else {
- sparse_cholesky = EigenSparseCholesky::Create(ordering_type);
+ sparse_cholesky = EigenSparseCholesky::Create(options.ordering_type);
}
break;
#else
@@ -77,9 +77,9 @@
case CX_SPARSE:
#ifndef CERES_NO_CXSPARSE
if (options.use_mixed_precision_solves) {
- sparse_cholesky = FloatCXSparseCholesky::Create(ordering_type);
+ sparse_cholesky = FloatCXSparseCholesky::Create(options.ordering_type);
} else {
- sparse_cholesky = CXSparseCholesky::Create(ordering_type);
+ sparse_cholesky = CXSparseCholesky::Create(options.ordering_type);
}
break;
#else
@@ -89,10 +89,11 @@
case ACCELERATE_SPARSE:
#ifndef CERES_NO_ACCELERATE_SPARSE
if (options.use_mixed_precision_solves) {
- sparse_cholesky = AppleAccelerateCholesky<float>::Create(ordering_type);
+ sparse_cholesky =
+ AppleAccelerateCholesky<float>::Create(options.ordering_type);
} else {
sparse_cholesky =
- AppleAccelerateCholesky<double>::Create(ordering_type);
+ AppleAccelerateCholesky<double>::Create(options.ordering_type);
}
break;
#else
diff --git a/internal/ceres/sparse_cholesky.h b/internal/ceres/sparse_cholesky.h
index f4f8996..feea7aa 100644
--- a/internal/ceres/sparse_cholesky.h
+++ b/internal/ceres/sparse_cholesky.h
@@ -63,7 +63,7 @@
// CompressedRowSparseMatrix lhs = ...;
// std::string message;
// CHECK_EQ(sparse_cholesky->Factorize(&lhs, &message),
-// LinearSolverTerminationType::SUCCESS);
+// LinearSolverTerminationType::SUCCESS);
// Vector rhs = ...;
// Vector solution = ...;
// CHECK_EQ(sparse_cholesky->Solve(rhs.data(), solution.data(), &message),
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc
index c9c8363..21b26a4 100644
--- a/internal/ceres/sparse_cholesky_test.cc
+++ b/internal/ceres/sparse_cholesky_test.cc
@@ -114,8 +114,7 @@
LinearSolver::Options sparse_cholesky_options;
sparse_cholesky_options.sparse_linear_algebra_library_type =
sparse_linear_algebra_library_type;
- sparse_cholesky_options.use_postordering =
- (ordering_type == OrderingType::AMD);
+ sparse_cholesky_options.ordering_type = ordering_type;
auto sparse_cholesky = SparseCholesky::Create(sparse_cholesky_options);
const CompressedRowSparseMatrix::StorageType storage_type =
sparse_cholesky->StorageType();
@@ -156,8 +155,7 @@
Param param = info.param;
std::stringstream ss;
ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
- << (::testing::get<1>(param) == OrderingType::AMD ? "AMD" : "NATURAL")
- << "_"
+ << ::testing::get<1>(param) << "_"
<< (::testing::get<2>(param) ? "UseBlockStructure" : "NoBlockStructure");
return ss.str();
}
@@ -197,8 +195,9 @@
SuiteSparseCholesky,
SparseCholeskyTest,
::testing::Combine(::testing::Values(SUITE_SPARSE),
- ::testing::Values(OrderingType::AMD,
- OrderingType::NATURAL),
+ ::testing::Values(OrderingType::NATURAL,
+ OrderingType::AMD,
+ OrderingType::NESDIS),
::testing::Values(true, false)),
ParamInfoToString);
#endif
diff --git a/internal/ceres/sparse_normal_cholesky_solver_test.cc b/internal/ceres/sparse_normal_cholesky_solver_test.cc
index 6b3316c..39d9d82 100644
--- a/internal/ceres/sparse_normal_cholesky_solver_test.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver_test.cc
@@ -111,7 +111,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = SUITE_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = false;
+ options.ordering_type = OrderingType::NATURAL;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -122,7 +122,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = SUITE_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = true;
+ options.ordering_type = OrderingType::AMD;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -135,7 +135,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = CX_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = false;
+ options.ordering_type = OrderingType::NATURAL;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -146,7 +146,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = CX_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = true;
+ options.ordering_type = OrderingType::AMD;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -159,7 +159,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = ACCELERATE_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = false;
+ options.ordering_type = OrderingType::NATURAL;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -170,7 +170,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = ACCELERATE_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = true;
+ options.ordering_type = OrderingType::AMD;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -183,7 +183,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = false;
+ options.ordering_type = OrderingType::NATURAL;
ContextImpl context;
options.context = &context;
TestSolver(options);
@@ -194,7 +194,7 @@
LinearSolver::Options options;
options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
options.type = SPARSE_NORMAL_CHOLESKY;
- options.use_postordering = true;
+ options.ordering_type = OrderingType::AMD;
ContextImpl context;
options.context = &context;
TestSolver(options);
diff --git a/internal/ceres/subset_preconditioner.cc b/internal/ceres/subset_preconditioner.cc
index c804274..a228545 100644
--- a/internal/ceres/subset_preconditioner.cc
+++ b/internal/ceres/subset_preconditioner.cc
@@ -51,7 +51,7 @@
LinearSolver::Options sparse_cholesky_options;
sparse_cholesky_options.sparse_linear_algebra_library_type =
options_.sparse_linear_algebra_library_type;
- sparse_cholesky_options.use_postordering = options_.use_postordering;
+ sparse_cholesky_options.ordering_type = options_.ordering_type;
sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options);
}
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index a7d28e6..82d20d4 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -43,6 +43,23 @@
#include "cholmod.h"
namespace ceres::internal {
+namespace {
+int OrderingTypeToCHOLMODEnum(OrderingType ordering_type) {
+ if (ordering_type == OrderingType::AMD) {
+ return CHOLMOD_AMD;
+ }
+ if (ordering_type == OrderingType::NESDIS) {
+ return CHOLMOD_NESDIS;
+ }
+
+ if (ordering_type == OrderingType::NATURAL) {
+ return CHOLMOD_NATURAL;
+ }
+ LOG(FATAL) << "Congratulations you have discovered a bug in Ceres Solver."
+ << "Please report it to the developers. " << ordering_type;
+ return -1;
+}
+} // namespace
using std::string;
using std::vector;
@@ -145,42 +162,38 @@
}
cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A,
+ OrderingType ordering_type,
string* message) {
- // Cholmod can try multiple re-ordering strategies to find a fill
- // reducing ordering. Here we just tell it use AMD with automatic
- // matrix dependence choice of supernodal versus simplicial
- // factorization.
cc_.nmethods = 1;
- cc_.method[0].ordering = CHOLMOD_AMD;
+ cc_.method[0].ordering = OrderingTypeToCHOLMODEnum(ordering_type);
cc_.supernodal = CHOLMOD_AUTO;
- cholmod_factor* factor = cholmod_analyze(A, &cc_);
- if (VLOG_IS_ON(2)) {
- cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
+ // TODO(sameeragarwal): The following modification exists entirely
+ // to preserve existing behaviour. It is not clear if this is
+ // actually the optimal thing to do. Once I understand the use of
+ // postordering better, I will come back and document this/remove it
+ // as the case maybe.
+ if (ordering_type == OrderingType::NATURAL) {
+ cc_.postorder = 0;
}
+ cholmod_factor* factor = cholmod_analyze(A, &cc_);
if (cc_.status != CHOLMOD_OK) {
*message =
StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
return nullptr;
}
-
CHECK(factor != nullptr);
- return factor;
-}
-cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(cholmod_sparse* A,
- const vector<int>& row_blocks,
- const vector<int>& col_blocks,
- string* message) {
- vector<int> ordering;
- if (!BlockAMDOrdering(A, row_blocks, col_blocks, &ordering)) {
- return nullptr;
+ if (VLOG_IS_ON(2)) {
+ cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
}
- return AnalyzeCholeskyWithUserOrdering(A, ordering, message);
}
-cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(
+return factor;
+}
+
+cholmod_factor* SuiteSparse::AnalyzeCholeskyWithGivenOrdering(
cholmod_sparse* A, const vector<int>& ordering, string* message) {
CHECK_EQ(ordering.size(), A->nrow);
@@ -189,9 +202,6 @@
cholmod_factor* factor =
cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), nullptr, 0, &cc_);
- if (VLOG_IS_ON(2)) {
- cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
- }
if (cc_.status != CHOLMOD_OK) {
*message =
StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
@@ -199,33 +209,26 @@
}
CHECK(factor != nullptr);
- return factor;
-}
-
-cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering(
- cholmod_sparse* A, string* message) {
- cc_.nmethods = 1;
- cc_.method[0].ordering = CHOLMOD_NATURAL;
- cc_.postorder = 0;
-
- cholmod_factor* factor = cholmod_analyze(A, &cc_);
if (VLOG_IS_ON(2)) {
cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
}
- if (cc_.status != CHOLMOD_OK) {
- *message =
- StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
- return nullptr;
- }
- CHECK(factor != nullptr);
return factor;
}
-bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A,
- const vector<int>& row_blocks,
- const vector<int>& col_blocks,
- vector<int>* ordering) {
+bool SuiteSparse::BlockOrdering(const cholmod_sparse* A,
+ OrderingType ordering_type,
+ const vector<int>& row_blocks,
+ const vector<int>& col_blocks,
+ vector<int>* ordering) {
+ if (ordering_type == OrderingType::NATURAL) {
+ ordering->resize(A->nrow);
+ for (int i = 0; i < A->nrow; ++i) {
+ (*ordering)[i] = i;
+ }
+ return true;
+ }
+
const int num_row_blocks = row_blocks.size();
const int num_col_blocks = col_blocks.size();
@@ -255,7 +258,7 @@
block_matrix.packed = 1;
vector<int> block_ordering(num_row_blocks);
- if (!cholmod_amd(&block_matrix, nullptr, 0, &block_ordering[0], &cc_)) {
+ if (!Ordering(&block_matrix, ordering_type, block_ordering.data())) {
return false;
}
@@ -263,6 +266,18 @@
return true;
}
+cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(cholmod_sparse* A,
+ OrderingType ordering_type,
+ const vector<int>& row_blocks,
+ const vector<int>& col_blocks,
+ string* message) {
+ vector<int> ordering;
+ if (!BlockOrdering(A, ordering_type, row_blocks, col_blocks, &ordering)) {
+ return nullptr;
+ }
+ return AnalyzeCholeskyWithGivenOrdering(A, ordering, message);
+}
+
LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A,
cholmod_factor* L,
string* message) {
@@ -333,13 +348,14 @@
return cholmod_solve(CHOLMOD_A, L, b, &cc_);
}
-bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
- int* ordering) {
- return cholmod_amd(matrix, nullptr, 0, ordering, &cc_);
-}
+bool SuiteSparse::Ordering(cholmod_sparse* matrix,
+ OrderingType ordering_type,
+ int* ordering) {
+ CHECK_NE(ordering_type, OrderingType::NATURAL);
+ if (ordering_type == OrderingType::AMD) {
+ return cholmod_amd(matrix, nullptr, 0, ordering, &cc_);
+ }
-bool SuiteSparse::NestedDissectionOrdering(cholmod_sparse* matrix,
- int* ordering) {
#ifdef CERES_NO_METIS
return false;
#else
@@ -355,6 +371,14 @@
return cholmod_camd(matrix, nullptr, 0, constraints, ordering, &cc_);
}
+bool SuiteSparse::IsNestedDissectionAvailable() {
+#ifdef CERES_NO_METIS
+ return false;
+#else
+ return true;
+#endif
+}
+
std::unique_ptr<SparseCholesky> SuiteSparseCholesky::Create(
const OrderingType ordering_type) {
return std::unique_ptr<SparseCholesky>(
@@ -379,23 +403,36 @@
cholmod_sparse cholmod_lhs = ss_.CreateSparseMatrixTransposeView(lhs);
+ // If a factorization does not exist, compute the symbolic
+ // factorization first.
+ //
+ // If the ordering type is NATURAL, then there is no fill reducing
+ // ordering to be computed, regardless of block structure, so we can
+ // just call the scalar version of symbolic factorization. For
+ // SuiteSparse this is the common case since we have already
+ // pre-ordered the columns of the Jacobian.
+ //
+ // Similarly regardless of ordering type, if there is no block
+ // structure in the matrix we call the scalar version of symbolic
+ // factorization.
if (factor_ == nullptr) {
- if (ordering_type_ == OrderingType::NATURAL) {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&cholmod_lhs, message);
+ if (ordering_type_ == OrderingType::NATURAL ||
+ (lhs->col_blocks().empty() || lhs->row_blocks().empty())) {
+ factor_ = ss_.AnalyzeCholesky(&cholmod_lhs, ordering_type_, message);
} else {
- if (!lhs->col_blocks().empty() && !(lhs->row_blocks().empty())) {
- factor_ = ss_.BlockAnalyzeCholesky(
- &cholmod_lhs, lhs->col_blocks(), lhs->row_blocks(), message);
- } else {
- factor_ = ss_.AnalyzeCholesky(&cholmod_lhs, message);
- }
- }
-
- if (factor_ == nullptr) {
- return LinearSolverTerminationType::FATAL_ERROR;
+ factor_ = ss_.BlockAnalyzeCholesky(&cholmod_lhs,
+ ordering_type_,
+ lhs->col_blocks(),
+ lhs->row_blocks(),
+ message);
}
}
+ if (factor_ == nullptr) {
+ return LinearSolverTerminationType::FATAL_ERROR;
+ }
+
+ // Compute and return the numeric factorization.
return ss_.Cholesky(&cholmod_lhs, factor_, message);
}
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index b8dc666..e0c485d 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -119,12 +119,11 @@
cholmod_sdmult(A, 0, alpha_, beta_, x, y, &cc_);
}
- // Find an ordering of A or AA' (if A is unsymmetric) that minimizes
- // the fill-in in the Cholesky factorization of the corresponding
- // matrix. This is done by using the AMD algorithm.
- //
- // Using this ordering, the symbolic Cholesky factorization of A (or
- // AA') is computed and returned.
+ // Compute a symbolic factorization for A or AA' (if A is
+ // unsymmetric). If ordering_type is NATURAL, then no fill reducing
+ // ordering is computed, otherwise depending on the value of
+ // ordering_type AMD or Nested Dissection is used to compute a fill
+ // reducing ordering before the symbolic factorization is computed.
//
// A is not modified, only the pattern of non-zeros of A is used,
// the actual numerical values in A are of no consequence.
@@ -132,9 +131,13 @@
// message contains an explanation of the failures if any.
//
// Caller owns the result.
- cholmod_factor* AnalyzeCholesky(cholmod_sparse* A, std::string* message);
+ cholmod_factor* AnalyzeCholesky(cholmod_sparse* A,
+ OrderingType ordering_type,
+ std::string* message);
+ // Block oriented version of AnalyzeCholesky.
cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A,
+ OrderingType ordering_type,
const std::vector<int>& row_blocks,
const std::vector<int>& col_blocks,
std::string* message);
@@ -150,20 +153,11 @@
// message contains an explanation of the failures if any.
//
// Caller owns the result.
- cholmod_factor* AnalyzeCholeskyWithUserOrdering(
+ cholmod_factor* AnalyzeCholeskyWithGivenOrdering(
cholmod_sparse* A,
const std::vector<int>& ordering,
std::string* message);
- // Perform a symbolic factorization of A without re-ordering A. No
- // postordering of the elimination tree is performed. This ensures
- // that the symbolic factor does not introduce an extra permutation
- // on the matrix. See the documentation for CHOLMOD for more details.
- //
- // message contains an explanation of the failures if any.
- cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A,
- std::string* message);
-
// Use the symbolic factorization in L, to find the numerical
// factorization for the matrix A or AA^T. Return true if
// successful, false otherwise. L contains the numeric factorization
@@ -183,50 +177,39 @@
cholmod_dense* b,
std::string* message);
+ // Find a fill reducing ordering. ordering is expected to be large
+ // enough to hold the ordering. ordering_type must be AMD or NESDIS.
+ bool Ordering(cholmod_sparse* matrix,
+ OrderingType ordering_type,
+ int* ordering);
+
+ // Find the block oriented fill reducing ordering of a matrix A,
+ // whose row and column blocks are given by row_blocks, and
+ // col_blocks respectively. The matrix may or may not be
+ // symmetric. The entries of col_blocks do not need to sum to the
+ // number of columns in A. If this is the case, only the first
+ // sum(col_blocks) are used to compute the ordering.
+ //
// By virtue of the modeling layer in Ceres being block oriented,
// all the matrices used by Ceres are also block oriented. When
// doing sparse direct factorization of these matrices the
- // fill-reducing ordering algorithms (in particular AMD) can either
- // be run on the block or the scalar form of these matrices. The two
- // SuiteSparse::AnalyzeCholesky methods allows the client to
- // compute the symbolic factorization of a matrix by either using
- // AMD on the matrix or a user provided ordering of the rows.
- //
- // But since the underlying matrices are block oriented, it is worth
- // running AMD on just the block structure of these matrices and then
- // lifting these block orderings to a full scalar ordering. This
- // preserves the block structure of the permuted matrix, and exposes
- // more of the super-nodal structure of the matrix to the numerical
- // factorization routines.
- //
- // Find the block oriented AMD ordering of a matrix A, whose row and
- // column blocks are given by row_blocks, and col_blocks
- // respectively. The matrix may or may not be symmetric. The entries
- // of col_blocks do not need to sum to the number of columns in
- // A. If this is the case, only the first sum(col_blocks) are used
- // to compute the ordering.
- bool BlockAMDOrdering(const cholmod_sparse* A,
- const std::vector<int>& row_blocks,
- const std::vector<int>& col_blocks,
- std::vector<int>* ordering);
-
- // Find a fill reducing approximate minimum degree
- // ordering. ordering is expected to be large enough to hold the
- // ordering.
- bool ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering);
-
- // Find a fill reducing ordering using nested dissection.
- bool NestedDissectionOrdering(cholmod_sparse* matrix, int* ordering);
+ // fill-reducing ordering algorithms can either be run on the block
+ // or the scalar form of these matrices. But since the underlying
+ // matrices are block oriented, it is worth running the fill
+ // reducing ordering on just the block structure of these matrices
+ // and then lifting these block orderings to a full scalar
+ // ordering. This preserves the block structure of the permuted
+ // matrix, and exposes more of the super-nodal structure of the
+ // matrix to the numerical factorization routines.
+ bool BlockOrdering(const cholmod_sparse* A,
+ OrderingType ordering_type,
+ const std::vector<int>& row_blocks,
+ const std::vector<int>& col_blocks,
+ std::vector<int>* ordering);
// Nested dissection is only available if SuiteSparse is compiled
// with Metis support.
- static bool IsNestedDissectionAvailable() {
-#ifdef CERES_NO_METIS
- return false;
-#else
- return true;
-#endif
- }
+ static bool IsNestedDissectionAvailable();
// Find a fill reducing approximate minimum degree
// ordering. constraints is an array which associates with each
@@ -300,6 +283,9 @@
class CERES_NO_EXPORT SuiteSparse {
public:
+ // Nested dissection is only available if SuiteSparse is compiled
+ // with Metis support.
+ static bool IsNestedDissectionAvailable() { return false; }
void Free(void* /*arg*/) {}
};
diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc
index daf6601..9467b5a 100644
--- a/internal/ceres/trust_region_preprocessor.cc
+++ b/internal/ceres/trust_region_preprocessor.cc
@@ -210,7 +210,10 @@
pp->linear_solver_options.max_num_refinement_iterations =
options.max_num_refinement_iterations;
pp->linear_solver_options.num_threads = options.num_threads;
- pp->linear_solver_options.use_postordering = options.use_postordering;
+ pp->linear_solver_options.ordering_type = OrderingType::NATURAL;
+ if (options.use_postordering) {
+ pp->linear_solver_options.ordering_type = OrderingType::AMD;
+ }
pp->linear_solver_options.context = pp->problem->context();
if (IsSchurType(pp->linear_solver_options.type)) {
@@ -236,7 +239,7 @@
// TODO(sameeragarwal): Implement the reordering of parameter
// blocks for CX_SPARSE.
if (options.sparse_linear_algebra_library_type == CX_SPARSE) {
- pp->linear_solver_options.use_postordering = true;
+ pp->linear_solver_options.ordering_type = OrderingType::AMD;
}
}
}
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index f04442d..bde2814 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -106,14 +106,7 @@
LinearSolver::Options sparse_cholesky_options;
sparse_cholesky_options.sparse_linear_algebra_library_type =
options_.sparse_linear_algebra_library_type;
-
- // The preconditioner's sparsity is not available in the
- // preprocessor, so the columns of the Jacobian have not been
- // reordered to minimize fill in when computing its sparse Cholesky
- // factorization. So we must tell the SparseCholesky object to
- // perform approximate minimum-degree reordering, which is done by
- // setting use_postordering to true.
- sparse_cholesky_options.use_postordering = true;
+ sparse_cholesky_options.ordering_type = options.ordering_type;
sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options);
const time_t init_time = time(nullptr);
