Enable mixed precision solves.
1. Add Solver::Options::use_mixed_precision_solves,
and Solver::Options::max_num_refinement_iterations.
2. Make SparseCholesky::Create return a unique_ptr.
3. SparseCholesky::Create now takes LinearSolver::Options
as an argument.
4. IterativeRefiner's constructor does not require num_cols
as an argument.
5. SparseNormalCholeskySolver now uses a separate rhs vector.
This basic implementation results in a 10% reduction in solver time
and 30% reduction in linear solver memory usage.
Change-Id: I6830f32cae2febf082d2733262eb2c9f0482b0ea
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 319e385..5619c52 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -120,6 +120,8 @@
"the pertubations.");
DEFINE_bool(line_search, false, "Use a line search instead of trust region "
"algorithm.");
+DEFINE_bool(mixed_precision_solves, false, "Use mixed precision solves.");
+DEFINE_int32(max_num_refinement_iterations, 0, "Iterative refinement iterations");
DEFINE_string(initial_ply, "", "Export the BAL file data as a PLY file.");
DEFINE_string(final_ply, "", "Export the refined BAL file data as a PLY "
"file.");
@@ -141,6 +143,8 @@
FLAGS_dense_linear_algebra_library,
&options->dense_linear_algebra_library_type));
options->use_explicit_schur_complement = FLAGS_explicit_schur_complement;
+ options->use_mixed_precision_solves = FLAGS_mixed_precision_solves;
+ options->max_num_refinement_iterations = FLAGS_max_num_refinement_iterations;
}
void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 592d329..9fc9ec5 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -497,9 +497,33 @@
// then it is probably best to keep this false, otherwise it will
// likely lead to worse performance.
- // This settings affects the SPARSE_NORMAL_CHOLESKY solver.
+ // This settings only affects the SPARSE_NORMAL_CHOLESKY solver.
bool dynamic_sparsity = false;
+ // TODO(sameeragarwal): Further expand the documentation for the
+ // following two options.
+
+ // NOTE1: EXPERIMETAL FEATURE, UNDER DEVELOPMENT, USE AT YOUR OWN RISK.
+ //
+ // If use_mixed_precision_solves is true, the Gauss-Newton matrix
+ // is computed in double precision, but its factorization is
+ // computed in single precision. This can result in significant
+ // time and memory savings at the cost of some accuracy in the
+ // Gauss-Newton step. Iterative refinement is used to recover some
+ // of this accuracy back.
+ //
+ // If use_mixed_precision_solves is true, we recommend setting
+ // max_num_refinement_iterations to 2-3.
+ //
+ // NOTE2: The following two options are currently only applicable
+ // if sparse_linear_algebra_library_type is EIGEN_SPARSE and
+ // linear_solver_type is SPARSE_NORMAL_CHOLESKY, or SPARSE_SCHUR.
+ bool use_mixed_precision_solves = false;
+
+ // Number steps of the iterative refinement process to run when
+ // computing the Gauss-Newton step.
+ int max_num_refinement_iterations = 0;
+
// Some non-linear least squares problems have additional
// structure in the way the parameter blocks interact that it is
// beneficial to modify the way the trust region step is computed.
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 41dbdff..8b57956 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -266,6 +266,8 @@
std::fill(values_.begin(), values_.end(), 0);
}
+// TODO(sameeragarwal): Make RightMultiply and LeftMultiply
+// block-aware for higher performance.
void CompressedRowSparseMatrix::RightMultiply(const double* x,
double* y) const {
CHECK_NOTNULL(x);
diff --git a/internal/ceres/cxsparse.cc b/internal/ceres/cxsparse.cc
index 11ff943..5a02877 100644
--- a/internal/ceres/cxsparse.cc
+++ b/internal/ceres/cxsparse.cc
@@ -196,8 +196,9 @@
void CXSparse::Free(csn* numeric_factor) { cs_di_nfree(numeric_factor); }
-CXSparseCholesky* CXSparseCholesky::Create(const OrderingType ordering_type) {
- return new CXSparseCholesky(ordering_type);
+std::unique_ptr<SparseCholesky> CXSparseCholesky::Create(
+ const OrderingType ordering_type) {
+ return std::unique_ptr<SparseCholesky>(new CXSparseCholesky(ordering_type));
}
CompressedRowSparseMatrix::StorageType CXSparseCholesky::StorageType() const {
diff --git a/internal/ceres/cxsparse.h b/internal/ceres/cxsparse.h
index 1ad79f9..1789afd 100644
--- a/internal/ceres/cxsparse.h
+++ b/internal/ceres/cxsparse.h
@@ -36,6 +36,7 @@
#ifndef CERES_NO_CXSPARSE
+#include <memory>
#include <string>
#include <vector>
@@ -140,7 +141,7 @@
class CXSparseCholesky : public SparseCholesky {
public:
// Factory
- static CXSparseCholesky* Create(const OrderingType ordering_type);
+ static std::unique_ptr<SparseCholesky> Create(OrderingType ordering_type);
// SparseCholesky interface.
virtual ~CXSparseCholesky();
diff --git a/internal/ceres/eigensparse.cc b/internal/ceres/eigensparse.cc
index 8ca1d17..425cd8b 100644
--- a/internal/ceres/eigensparse.cc
+++ b/internal/ceres/eigensparse.cc
@@ -41,6 +41,8 @@
namespace ceres {
namespace internal {
+// TODO(sameeragarwal): Use enable_if to clean up the implementations
+// for when Scalar == double.
template <typename Solver>
class EigenSparseCholeskyTemplate : public SparseCholesky {
public:
@@ -78,23 +80,22 @@
return LINEAR_SOLVER_SUCCESS;
}
- virtual LinearSolverTerminationType Solve(const double* rhs_ptr,
- double* solution_ptr,
- std::string* message) {
+ LinearSolverTerminationType Solve(const double* rhs_ptr,
+ double* solution_ptr,
+ std::string* message) {
CHECK(analyzed_) << "Solve called without a call to Factorize first.";
- ConstVectorRef rhs(rhs_ptr, solver_.cols());
- VectorRef solution(solution_ptr, solver_.cols());
+ scalar_rhs_ = ConstVectorRef(rhs_ptr, solver_.cols())
+ .template cast<typename Solver::Scalar>();
// The two casts are needed if the Scalar in this class is not
// double. For code simplicitly we are going to assume that Eigen
// is smart enough to figure out that casting a double Vector to a
// double Vector is a straight copy. If this turns into a
// performance bottleneck (unlikely), we can revisit this.
- solution =
- solver_
- .solve(rhs.template cast<typename Solver::Scalar>())
- .template cast<double>();
+ scalar_solution_ = solver_.solve(scalar_rhs_);
+ VectorRef(solution_ptr, solver_.cols()) =
+ scalar_solution_.template cast<double>();
if (solver_.info() != Eigen::Success) {
*message = "Eigen failure. Unable to do triangular solve.";
@@ -131,12 +132,16 @@
}
private:
- Eigen::Matrix<typename Solver::Scalar, Eigen::Dynamic, 1> values_;
+ Eigen::Matrix<typename Solver::Scalar, Eigen::Dynamic, 1> values_,
+ scalar_rhs_, scalar_solution_;
bool analyzed_;
Solver solver_;
};
-SparseCholesky* EigenSparseCholesky::Create(const OrderingType ordering_type) {
+std::unique_ptr<SparseCholesky> EigenSparseCholesky::Create(
+ const OrderingType ordering_type) {
+ std::unique_ptr<SparseCholesky> sparse_cholesky;
+
// The preprocessor gymnastics here are dealing with the fact that
// before version 3.2.2, Eigen did not support a third template
// parameter to specify the ordering and it always defaults to AMD.
@@ -150,44 +155,48 @@
Eigen::NaturalOrdering<int>>
WithNaturalOrdering;
if (ordering_type == AMD) {
- return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
+ sparse_cholesky.reset(new EigenSparseCholeskyTemplate<WithAMDOrdering>());
} else {
- return new EigenSparseCholeskyTemplate<WithNaturalOrdering>();
+ sparse_cholesky.reset(
+ new EigenSparseCholeskyTemplate<WithNaturalOrdering>());
}
#else
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper>
WithAMDOrdering;
- return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
+ sparse_cholesky.reset(new EigenSparseCholeskyTemplate<WithAMDOrdering>());
#endif
+ return sparse_cholesky;
}
EigenSparseCholesky::~EigenSparseCholesky() {}
-SparseCholesky* FloatEigenSparseCholesky::Create(
+std::unique_ptr<SparseCholesky> FloatEigenSparseCholesky::Create(
const OrderingType ordering_type) {
+ std::unique_ptr<SparseCholesky> sparse_cholesky;
// The preprocessor gymnastics here are dealing with the fact that
// before version 3.2.2, Eigen did not support a third template
// parameter to specify the ordering and it always defaults to AMD.
#if EIGEN_VERSION_AT_LEAST(3, 2, 2)
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>,
Eigen::Upper,
- Eigen::AMDOrdering<int> >
+ Eigen::AMDOrdering<int>>
WithAMDOrdering;
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>,
Eigen::Upper,
- Eigen::NaturalOrdering<int> >
+ Eigen::NaturalOrdering<int>>
WithNaturalOrdering;
if (ordering_type == AMD) {
- LOG(FATAL) << "We should not be doing this";
- return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
+ sparse_cholesky.reset(new EigenSparseCholeskyTemplate<WithAMDOrdering>());
} else {
- return new EigenSparseCholeskyTemplate<WithNaturalOrdering>();
+ sparse_cholesky.reset(
+ new EigenSparseCholeskyTemplate<WithNaturalOrdering>());
}
#else
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>, Eigen::Upper>
WithAMDOrdering;
- return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
+ sparse_cholesky.reset(new EigenSparseCholeskyTemplate<WithAMDOrdering>());
#endif
+ return sparse_cholesky;
}
FloatEigenSparseCholesky::~FloatEigenSparseCholesky() {}
diff --git a/internal/ceres/eigensparse.h b/internal/ceres/eigensparse.h
index ce8ef34..2e6c6f0 100644
--- a/internal/ceres/eigensparse.h
+++ b/internal/ceres/eigensparse.h
@@ -38,7 +38,9 @@
#ifdef CERES_USE_EIGEN_SPARSE
+#include <memory>
#include <string>
+
#include "Eigen/SparseCore"
#include "ceres/linear_solver.h"
#include "ceres/sparse_cholesky.h"
@@ -49,7 +51,8 @@
class EigenSparseCholesky : public SparseCholesky {
public:
// Factory
- static SparseCholesky* Create(const OrderingType ordering_type);
+ static std::unique_ptr<SparseCholesky> Create(
+ const OrderingType ordering_type);
// SparseCholesky interface.
virtual ~EigenSparseCholesky();
@@ -66,7 +69,8 @@
class FloatEigenSparseCholesky : public SparseCholesky {
public:
// Factory
- static SparseCholesky* Create(const OrderingType ordering_type);
+ static std::unique_ptr<SparseCholesky> Create(
+ const OrderingType ordering_type);
// SparseCholesky interface.
virtual ~FloatEigenSparseCholesky();
diff --git a/internal/ceres/iterative_refiner.cc b/internal/ceres/iterative_refiner.cc
index fff343d..b73c0b0 100644
--- a/internal/ceres/iterative_refiner.cc
+++ b/internal/ceres/iterative_refiner.cc
@@ -38,25 +38,28 @@
namespace ceres {
namespace internal {
-IterativeRefiner::IterativeRefiner(const int num_cols,
- const int max_num_iterations)
- : num_cols_(num_cols),
- max_num_iterations_(max_num_iterations),
- residual_(num_cols),
- correction_(num_cols),
- lhs_x_solution_(num_cols) {}
+IterativeRefiner::IterativeRefiner(const int max_num_iterations)
+ : max_num_iterations_(max_num_iterations) {}
IterativeRefiner::~IterativeRefiner() {}
+void IterativeRefiner::Allocate(int num_cols) {
+ residual_.resize(num_cols);
+ correction_.resize(num_cols);
+ lhs_x_solution_.resize(num_cols);
+}
+
IterativeRefiner::Summary IterativeRefiner::Refine(
const SparseMatrix& lhs,
const double* rhs_ptr,
SparseCholesky* sparse_cholesky,
double* solution_ptr) {
Summary summary;
+ const int num_cols = lhs.num_cols();
+ Allocate(num_cols);
- ConstVectorRef rhs(rhs_ptr, num_cols_);
- VectorRef solution(solution_ptr, num_cols_);
+ ConstVectorRef rhs(rhs_ptr, num_cols);
+ VectorRef solution(solution_ptr, num_cols);
summary.lhs_max_norm = ConstVectorRef(lhs.values(), lhs.num_nonzeros())
.lpNorm<Eigen::Infinity>();
diff --git a/internal/ceres/iterative_refiner.h b/internal/ceres/iterative_refiner.h
index d1efb5c..9021e81 100644
--- a/internal/ceres/iterative_refiner.h
+++ b/internal/ceres/iterative_refiner.h
@@ -77,12 +77,9 @@
double residual_max_norm = -1;
};
- // num_cols is the number of rows & columns in the linear system
- // being solved.
- //
// max_num_iterations is the maximum number of refinement iterations
// to perform.
- IterativeRefiner(int num_cols, int max_num_iterations);
+ IterativeRefiner(int max_num_iterations);
// Needed for mocking.
virtual ~IterativeRefiner();
@@ -97,13 +94,18 @@
// to lhs * x = rhs. It can be zero.
//
// This method is virtual to facilitate mocking.
+ //
+ // TODO(sameeragarwal): Consider dropping the Summary object, and
+ // simplifying the internal implementation to improve efficiency,
+ // since we do not seem to be using the output at all.
virtual Summary Refine(const SparseMatrix& lhs,
const double* rhs,
SparseCholesky* sparse_cholesky,
double* solution);
private:
- int num_cols_;
+ void Allocate(int num_cols);
+
int max_num_iterations_;
Vector residual_;
Vector correction_;
diff --git a/internal/ceres/iterative_refiner_test.cc b/internal/ceres/iterative_refiner_test.cc
index a3bcbf1..0c9ef69 100644
--- a/internal/ceres/iterative_refiner_test.cc
+++ b/internal/ceres/iterative_refiner_test.cc
@@ -147,7 +147,7 @@
ExactSolutionWithExactFactorizationReturnsInZeroIterations) {
FakeSparseMatrix lhs(lhs_);
FakeSparseCholesky<double> sparse_cholesky(lhs_);
- IterativeRefiner refiner(num_cols_, max_num_iterations_);
+ IterativeRefiner refiner(max_num_iterations_);
Vector refined_solution = solution_;
auto summary = refiner.Refine(
lhs, rhs_.data(), &sparse_cholesky, refined_solution.data());
@@ -161,7 +161,7 @@
RandomSolutionWithExactFactorizationReturnsInOneIteration) {
FakeSparseMatrix lhs(lhs_);
FakeSparseCholesky<double> sparse_cholesky(lhs_);
- IterativeRefiner refiner(num_cols_, max_num_iterations_);
+ IterativeRefiner refiner(max_num_iterations_);
Vector refined_solution(num_cols_);
refined_solution.setRandom();
auto summary = refiner.Refine(
@@ -178,7 +178,7 @@
// Use a single precision Cholesky factorization of the double
// precision matrix. This will give us an approximate factorization.
FakeSparseCholesky<float> sparse_cholesky(lhs_);
- IterativeRefiner refiner(num_cols_, max_num_iterations_);
+ IterativeRefiner refiner(max_num_iterations_);
Vector refined_solution(num_cols_);
refined_solution.setRandom();
auto summary = refiner.Refine(
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index a7c61bc..51897fc 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -160,6 +160,8 @@
int e_block_size = Eigen::Dynamic;
int f_block_size = Eigen::Dynamic;
+ bool use_mixed_precision_solves = false;
+ int max_num_refinement_iterations = 0;
ContextImpl* context = nullptr;
};
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 7934e36..e051b98 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -227,9 +227,7 @@
const LinearSolver::Options& options)
: SchurComplementSolver(options) {
if (options.type != ITERATIVE_SCHUR) {
- sparse_cholesky_.reset(
- SparseCholesky::Create(options.sparse_linear_algebra_library_type,
- options.use_postordering ? AMD : NATURAL));
+ sparse_cholesky_ = SparseCholesky::Create(options);
}
}
diff --git a/internal/ceres/sparse_cholesky.cc b/internal/ceres/sparse_cholesky.cc
index 3c4e97f..ef6d762 100644
--- a/internal/ceres/sparse_cholesky.cc
+++ b/internal/ceres/sparse_cholesky.cc
@@ -38,41 +38,53 @@
namespace ceres {
namespace internal {
-SparseCholesky* SparseCholesky::Create(
- SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
- OrderingType ordering_type) {
- switch (sparse_linear_algebra_library_type) {
+std::unique_ptr<SparseCholesky> SparseCholesky::Create(
+ const LinearSolver::Options& options) {
+ const OrderingType ordering_type = options.use_postordering ? AMD : NATURAL;
+ std::unique_ptr<SparseCholesky> sparse_cholesky;
+
+ switch (options.sparse_linear_algebra_library_type) {
case SUITE_SPARSE:
#ifndef CERES_NO_SUITESPARSE
- return SuiteSparseCholesky::Create(ordering_type);
+ sparse_cholesky = SuiteSparseCholesky::Create(ordering_type);
+ break;
#else
LOG(FATAL) << "Ceres was compiled without support for SuiteSparse.";
- return NULL;
#endif
case EIGEN_SPARSE:
#ifdef CERES_USE_EIGEN_SPARSE
- return EigenSparseCholesky::Create(ordering_type);
+ if (options.use_mixed_precision_solves) {
+ sparse_cholesky = FloatEigenSparseCholesky::Create(ordering_type);
+ } else {
+ sparse_cholesky = EigenSparseCholesky::Create(ordering_type);
+ }
+ break;
#else
LOG(FATAL) << "Ceres was compiled without support for "
<< "Eigen's sparse Cholesky factorization routines.";
- return NULL;
#endif
case CX_SPARSE:
#ifndef CERES_NO_CXSPARSE
- return CXSparseCholesky::Create(ordering_type);
+ sparse_cholesky = CXSparseCholesky::Create(ordering_type);
#else
LOG(FATAL) << "Ceres was compiled without support for CXSparse.";
- return NULL;
#endif
-
+ break;
default:
LOG(FATAL) << "Unknown sparse linear algebra library type : "
<< SparseLinearAlgebraLibraryTypeToString(
- sparse_linear_algebra_library_type);
+ options.sparse_linear_algebra_library_type);
}
- return NULL;
+
+ if (options.max_num_refinement_iterations > 0) {
+ std::unique_ptr<IterativeRefiner> refiner(
+ new IterativeRefiner(options.max_num_refinement_iterations));
+ sparse_cholesky = std::unique_ptr<SparseCholesky>(new RefinedSparseCholesky(
+ std::move(sparse_cholesky), std::move(refiner)));
+ }
+ return sparse_cholesky;
}
SparseCholesky::~SparseCholesky() {}
@@ -105,7 +117,8 @@
RefinedSparseCholesky::~RefinedSparseCholesky() {}
-CompressedRowSparseMatrix::StorageType RefinedSparseCholesky::StorageType() const {
+CompressedRowSparseMatrix::StorageType RefinedSparseCholesky::StorageType()
+ const {
return sparse_cholesky_->StorageType();
}
diff --git a/internal/ceres/sparse_cholesky.h b/internal/ceres/sparse_cholesky.h
index c0e3e86..bbe4237 100644
--- a/internal/ceres/sparse_cholesky.h
+++ b/internal/ceres/sparse_cholesky.h
@@ -66,14 +66,8 @@
class SparseCholesky {
public:
- // Factory which returns an instance of SparseCholesky for the given
- // sparse linear algebra library and fill reducing ordering
- // strategy.
- //
- // Caller owns the result.
- static SparseCholesky* Create(
- SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
- OrderingType ordering_type);
+ static std::unique_ptr<SparseCholesky> Create(
+ const LinearSolver::Options& options);
virtual ~SparseCholesky();
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc
index f0cd729..1c79bc9 100644
--- a/internal/ceres/sparse_cholesky_test.cc
+++ b/internal/ceres/sparse_cholesky_test.cc
@@ -109,8 +109,12 @@
const int min_block_size,
const int max_block_size,
const double block_density) {
- std::unique_ptr<SparseCholesky> sparse_cholesky(SparseCholesky::Create(
- sparse_linear_algebra_library_type, ordering_type));
+ 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 == AMD);
+ std::unique_ptr<SparseCholesky> sparse_cholesky = SparseCholesky::Create(
+ sparse_cholesky_options);
const CompressedRowSparseMatrix::StorageType storage_type =
sparse_cholesky->StorageType();
@@ -230,7 +234,7 @@
class MockIterativeRefiner : public IterativeRefiner {
public:
- MockIterativeRefiner() : IterativeRefiner(1, 1) {}
+ MockIterativeRefiner() : IterativeRefiner(1) {}
MOCK_METHOD4(Refine,
Summary(const SparseMatrix& lhs,
const double* rhs,
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 0572870..0f2e589 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -38,6 +38,7 @@
#include "ceres/block_sparse_matrix.h"
#include "ceres/inner_product_computer.h"
#include "ceres/internal/eigen.h"
+#include "ceres/iterative_refiner.h"
#include "ceres/linear_solver.h"
#include "ceres/sparse_cholesky.h"
#include "ceres/triplet_sparse_matrix.h"
@@ -50,9 +51,7 @@
SparseNormalCholeskySolver::SparseNormalCholeskySolver(
const LinearSolver::Options& options)
: options_(options) {
- sparse_cholesky_.reset(
- SparseCholesky::Create(options_.sparse_linear_algebra_library_type,
- options_.use_postordering ? AMD : NATURAL));
+ sparse_cholesky_ = SparseCholesky::Create(options);
}
SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {}
@@ -69,8 +68,11 @@
summary.message = "Success.";
const int num_cols = A->num_cols();
- VectorRef(x, num_cols).setZero();
- A->LeftMultiply(b, x);
+ VectorRef xref(x, num_cols);
+ xref.setZero();
+ rhs_.resize(num_cols);
+ rhs_.setZero();
+ A->LeftMultiply(b, rhs_.data());
event_logger.AddEvent("Compute RHS");
if (per_solve_options.D != NULL) {
@@ -95,14 +97,16 @@
inner_product_computer_->Compute();
event_logger.AddEvent("InnerProductComputer::Compute");
- // TODO(sameeragarwal):
-
if (per_solve_options.D != NULL) {
A->DeleteRowBlocks(A->block_structure()->cols.size());
}
+
summary.termination_type = sparse_cholesky_->FactorAndSolve(
- inner_product_computer_->mutable_result(), x, x, &summary.message);
- event_logger.AddEvent("Factor & Solve");
+ inner_product_computer_->mutable_result(),
+ rhs_.data(),
+ x,
+ &summary.message);
+ event_logger.AddEvent("SparseCholesky::FactorAndSolve");
return summary;
}
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index 537a2c9..59c7c0b 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -63,6 +63,7 @@
double* x);
const LinearSolver::Options options_;
+ Vector rhs_;
std::unique_ptr<SparseCholesky> sparse_cholesky_;
std::unique_ptr<InnerProductComputer> inner_product_computer_;
CERES_DISALLOW_COPY_AND_ASSIGN(SparseNormalCholeskySolver);
diff --git a/internal/ceres/subset_preconditioner.cc b/internal/ceres/subset_preconditioner.cc
index c7f5ef6..de6e9fa 100644
--- a/internal/ceres/subset_preconditioner.cc
+++ b/internal/ceres/subset_preconditioner.cc
@@ -44,10 +44,13 @@
SubsetPreconditioner::SubsetPreconditioner(
const Preconditioner::Options& options, const BlockSparseMatrix& A)
: options_(options), num_cols_(A.num_cols()) {
- sparse_cholesky_.reset(
- SparseCholesky::Create(options_.sparse_linear_algebra_library_type,
- options_.use_postordering ? AMD : NATURAL));
CHECK_GE(options_.subset_preconditioner_start_row_block, 0);
+ 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_ = SparseCholesky::Create(sparse_cholesky_options);
}
SubsetPreconditioner::~SubsetPreconditioner() {}
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index cc07a01..f10f57e 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -348,9 +348,9 @@
#endif
}
-SuiteSparseCholesky* SuiteSparseCholesky::Create(
+std::unique_ptr<SparseCholesky> SuiteSparseCholesky::Create(
const OrderingType ordering_type) {
- return new SuiteSparseCholesky(ordering_type);
+ return std::unique_ptr<SparseCholesky>(new SuiteSparseCholesky(ordering_type));
}
SuiteSparseCholesky::SuiteSparseCholesky(const OrderingType ordering_type)
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 91cd24b..6fef0fd 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -285,7 +285,8 @@
class SuiteSparseCholesky : public SparseCholesky {
public:
- static SuiteSparseCholesky* Create(const OrderingType ordering_type);
+ static std::unique_ptr<SparseCholesky> Create(
+ OrderingType ordering_type);
// SparseCholesky interface.
virtual ~SuiteSparseCholesky();
diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc
index 96f597c..cca2cf7 100644
--- a/internal/ceres/trust_region_preprocessor.cc
+++ b/internal/ceres/trust_region_preprocessor.cc
@@ -193,6 +193,10 @@
pp->linear_solver_options.use_explicit_schur_complement =
options.use_explicit_schur_complement;
pp->linear_solver_options.dynamic_sparsity = options.dynamic_sparsity;
+ pp->linear_solver_options.use_mixed_precision_solves =
+ options.use_mixed_precision_solves;
+ 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.context = pp->problem->context();
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index 36247b7..7cd9bd4 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -105,8 +105,18 @@
InitEliminator(bs);
const time_t eliminator_time = time(NULL);
- sparse_cholesky_.reset(
- SparseCholesky::Create(options_.sparse_linear_algebra_library_type, AMD));
+ 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 approximiate minimum-degree reordering, which is done by
+ // setting use_postordering to true.
+ sparse_cholesky_options.use_postordering = true;
+ sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options);
const time_t init_time = time(NULL);
VLOG(2) << "init time: " << init_time - start_time
