saberbud/ceres-solver
1
0
Fork 0
mirror of https://github.com/ceres-solver/ceres-solver.git synced 2026-08-29 08:34:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add mixed precision solves for SUITE_SPARSE

Browse Source 
Starting with SuiteSparse version 7.4.0 CHOLMOD has support for single
precision matrices. This allows us to have single precision and mixed
precision solves when using the SUITE_SPARSE backend.

This CL also fixes sparse_cholesky_test which was completely broken for
single precision testing.

Sample performance on my Mac.
/usr/bin/time -l ./bin/bundle_adjuster --input=../../Downloads/problem-3068-310854-pre.txt
<SNIP>

Cost:
Initial                          9.099334e+07
Final                            4.161838e+06
Change                           8.683150e+07

Minimizer iterations                        6
Successful steps                            4
Unsuccessful steps                          2

Time (in seconds):
Preprocessor                         2.528222

  Residual only evaluation           0.142804 (5)
  Jacobian & residual evaluation     0.424014 (4)
  Linear solver                     54.083396 (5)
Minimizer                           54.895752

Postprocessor                        0.024564
Total                               57.448539

Termination:                   NO_CONVERGENCE (Maximum number of iterations reached. Number of iterations: 5.)

       59.04 real       341.24 user         5.49 sys
          5776375808  maximum resident set size
<SNIP>
        616329634071  instructions retired
        929475980510  cycles elapsed
          5375034560  peak memory footprint

/usr/bin/time -l ./bin/bundle_adjuster --input=../../Downloads/problem-3068-310854-pre.txt  -mixed_precision_solves
<SNIP>

Cost:
Initial                          9.099334e+07
Final                            4.148930e+06
Change                           8.684441e+07

Minimizer iterations                        6
Successful steps                            4
Unsuccessful steps                          2

Time (in seconds):
Preprocessor                         2.580217

  Residual only evaluation           0.144098 (5)
  Jacobian & residual evaluation     0.396723 (4)
  Linear solver                     23.636074 (5)
Minimizer                           24.427163

Postprocessor                        0.023790
Total                               27.031170

Termination:                   NO_CONVERGENCE (Maximum number of iterations reached. Number of iterations: 5.)

       28.58 real       128.53 user         2.37 sys
          4818386944  maximum resident set size
<SNIP>
        395186936091  instructions retired
        368802808856  cycles elapsed
          4327029824  peak memory footprint

Change-Id: I1f137b0dd12da8da7f9ced338dd8f20f4bbdf99d
This commit is contained in:
Sameer Agarwal
2024-08-24 11:18:24 -07:00
parent 8f85014dcf
commit 80fce72bfd
13 changed files with 276 additions and 207 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+5
View File
@@ -374,6 +374,11 @@ if (SUITESPARSE)
message("-- Found SuiteSparse ${SuiteSparse_VERSION}, "
"building with SuiteSparse.")
if (SuiteSparse_VERSION VERSION_LESS 7.4.0)
list(APPEND CERES_COMPILE_OPTIONS CERES_NO_CHOLMOD_FLOAT)
else()
message("-- CHOLMOD supports single precision factorization")
endif()
if (SuiteSparse_NO_CMAKE OR NOT SuiteSparse_DIR)
install(FILES ${Ceres_SOURCE_DIR}/cmake/FindSuiteSparse.cmake
${Ceres_SOURCE_DIR}/cmake/FindMETIS.cmake
cmake/config.h.in
+5
View File
@@ -77,6 +77,11 @@
// If defined, Ceres was compiled with a version of SuiteSparse/CHOLMOD without
// the Partition module (requires METIS).
@CERES_NO_CHOLMOD_PARTITION@
// If defined Ceres was compiled with a version of SuiteSparse that does not
// support single precision solves in CHOLMOD.
@CERES_NO_CHOLMOD_FLOAT@
// If defined Ceres was compiled without support for METIS via Eigen.
@CERES_NO_EIGEN_METIS@
docs/source/nnls_solving.rst
+5 -13
View File
@@ -905,20 +905,12 @@ iterations of iterative refinement are controlled by
:member:`Solver::Options::max_num_refinement_iterations`. The default
value of this parameter is zero, which means if
:member:`Solver::Options::use_mixed_precision_solves` is ``true``,
then no iterative refinement is performed. Usually 2-3 refinement
iterations are enough.
then no iterative refinement is performed. Usually 1-3 refinement
iterations are enough, depending upon the conditioning of your
problem.
Mixed precision solves are available in the following linear solver
configurations:
1. ``DENSE_NORMAL_CHOLESKY`` + ``EIGEN``/ ``LAPACK`` / ``CUDA``.
2. ``DENSE_SCHUR`` + ``EIGEN``/ ``LAPACK`` / ``CUDA``.
3. ``SPARSE_NORMAL_CHOLESKY`` + ``EIGEN_SPARSE`` / ``ACCELERATE_SPARSE``
4. ``SPARSE_SCHUR`` + ``EIGEN_SPARSE`` / ``ACCELERATE_SPARSE``
Mixed precision solves area not available when using ``SUITE_SPARSE``
as the sparse linear algebra backend because SuiteSparse/CHOLMOD does
not support single precision solves.
If :member:`Solver::Options::max_num_refinement_iterations = 0`, then
the Gauss-Newton step is computed in single precision.
.. _section-preconditioner:
include/ceres/solver.h
+12 -11
View File
@@ -569,23 +569,24 @@ class CERES_EXPORT Solver {
// This settings only affects the SPARSE_NORMAL_CHOLESKY solver.
bool dynamic_sparsity = false;
// 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, 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.
// max_num_refinement_iterations to 1-3, depending on your problem. If
// max_num_refinement_iterations = 0, then the Gauss-Newton step is computed
// in single precision.
//
// This options is available when linear solver uses sparse or dense
// cholesky factorization, except when sparse_linear_algebra_library_type =
// SUITE_SPARSE.
// cholesky factorization.
bool use_mixed_precision_solves = false;
// Number steps of the iterative refinement process to run when
// computing the Gauss-Newton step.
// Number steps of the iterative refinement process to run when computing
// the Gauss-Newton step. This is most useful when used in conjunction with
// use_mixed_precision = true.
int max_num_refinement_iterations = 0;
// Minimum number of iterations for which the linear solver should
internal/ceres/CMakeLists.txt
-1
View File
@@ -184,7 +184,6 @@ add_library(ceres_internal OBJECT
fake_bundle_adjustment_jacobian.cc
file.cc
first_order_function.cc
float_suitesparse.cc
function_sample.cc
gradient_checking_cost_function.cc
gradient_problem_solver.cc
internal/ceres/float_suitesparse.cc
-49
View File
@@ -1,49 +0,0 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2023 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/float_suitesparse.h"
#include <memory>
#include "absl/log/log.h"
#if !defined(CERES_NO_SUITESPARSE)
namespace ceres::internal {
std::unique_ptr<SparseCholesky> FloatSuiteSparseCholesky::Create(
OrderingType ordering_type) {
LOG(FATAL) << "FloatSuiteSparseCholesky is not available.";
return {};
}
} // namespace ceres::internal
#endif // !defined(CERES_NO_SUITESPARSE)
internal/ceres/float_suitesparse.h
-59
View File
@@ -1,59 +0,0 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2023 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_FLOAT_SUITESPARSE_H_
#define CERES_INTERNAL_FLOAT_SUITESPARSE_H_
// This include must come before any #ifndef check on Ceres compile options.
// clang-format off
#include "ceres/internal/config.h"
// clang-format on
#include <memory>
#include "ceres/internal/export.h"
#include "ceres/sparse_cholesky.h"
#if !defined(CERES_NO_SUITESPARSE)
namespace ceres::internal {
// Fake implementation of a single precision Sparse Cholesky using
// SuiteSparse.
class CERES_NO_EXPORT FloatSuiteSparseCholesky : public SparseCholesky {
public:
static std::unique_ptr<SparseCholesky> Create(OrderingType ordering_type);
};
} // namespace ceres::internal
#endif // !defined(CERES_NO_SUITESPARSE)
#endif // CERES_INTERNAL_FLOAT_SUITESPARSE_H_
internal/ceres/solver.cc
+6 -1
View File
@@ -195,11 +195,16 @@ bool OptionsAreValidForSparseCholeskyBasedSolver(const Solver::Options& options,
return false;
}
#ifdef CERES_NO_CHOLMOD_FLOAT
if (options.use_mixed_precision_solves &&
options.sparse_linear_algebra_library_type == SUITE_SPARSE) {
*error = absl::StrFormat(kMixedFormat, solver_name, library_name);
*error =
"This version of SuiteSparse does not support single precision "
"Cholesky factorization. So use_mixed_precision_solves is not "
"supported with SUITE_SPARSE";
return false;
}
#endif
if (options.dynamic_sparsity &&
options.sparse_linear_algebra_library_type == ACCELERATE_SPARSE) {
internal/ceres/solver_test.cc
+24
View File
@@ -511,7 +511,11 @@ TEST(Solver, SparseNormalCholeskyOptionsSuiteSparse) {
options.sparse_linear_algebra_library_type)) {
options.use_mixed_precision_solves = true;
options.dynamic_sparsity = false;
#ifdef CERES_NO_CHOLMOD_FLOAT
EXPECT_FALSE(options.IsValid(&message));
#else
EXPECT_TRUE(options.IsValid(&message));
#endif
options.use_mixed_precision_solves = false;
options.dynamic_sparsity = true;
@@ -519,7 +523,11 @@ TEST(Solver, SparseNormalCholeskyOptionsSuiteSparse) {
options.use_mixed_precision_solves = true;
options.dynamic_sparsity = true;
#ifdef CERES_NO_CHOLMOD_FLOAT
EXPECT_FALSE(options.IsValid(&message));
#else
EXPECT_TRUE(options.IsValid(&message));
#endif
}
#ifndef CERES_NO_CHOLMOD_PARTITION
@@ -532,7 +540,11 @@ TEST(Solver, SparseNormalCholeskyOptionsSuiteSparse) {
options.use_mixed_precision_solves = true;
options.dynamic_sparsity = false;
#ifdef CERES_NO_CHOLMOD_FLOAT
EXPECT_FALSE(options.IsValid(&message));
#else
EXPECT_TRUE(options.IsValid(&message));
#endif
options.use_mixed_precision_solves = false;
options.dynamic_sparsity = true;
@@ -540,7 +552,11 @@ TEST(Solver, SparseNormalCholeskyOptionsSuiteSparse) {
options.use_mixed_precision_solves = true;
options.dynamic_sparsity = true;
#ifdef CERES_NO_CHOLMOD_FLOAT
EXPECT_FALSE(options.IsValid(&message));
#else
EXPECT_TRUE(options.IsValid(&message));
#endif
}
#else
options.linear_solver_ordering_type = NESDIS;
@@ -729,7 +745,11 @@ TEST(Solver, SparseSchurOptionsSuiteSparse) {
options.sparse_linear_algebra_library_type)) {
options.use_mixed_precision_solves = true;
options.dynamic_sparsity = false;
#ifdef CERES_NO_CHOLMOD_FLOAT
EXPECT_FALSE(options.IsValid(&message));
#else
EXPECT_TRUE(options.IsValid(&message));
#endif
options.use_mixed_precision_solves = false;
options.dynamic_sparsity = true;
@@ -750,7 +770,11 @@ TEST(Solver, SparseSchurOptionsSuiteSparse) {
options.use_mixed_precision_solves = true;
options.dynamic_sparsity = false;
#ifdef CERES_NO_CHOLMOD_FLOAT
EXPECT_FALSE(options.IsValid(&message));
#else
EXPECT_TRUE(options.IsValid(&message));
#endif
options.use_mixed_precision_solves = false;
options.dynamic_sparsity = true;
internal/ceres/sparse_cholesky.cc
-1
View File
@@ -38,7 +38,6 @@
#include "ceres/accelerate_sparse.h"
#include "ceres/cuda_sparse_cholesky.h"
#include "ceres/eigensparse.h"
#include "ceres/float_suitesparse.h"
#include "ceres/iterative_refiner.h"
#include "ceres/suitesparse.h"
internal/ceres/sparse_cholesky_test.cc
+45 -69
View File
@@ -113,6 +113,7 @@ bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
void SparseCholeskySolverUnitTest(
const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const bool use_single_precision,
const OrderingType ordering_type,
const bool use_block_structure,
const int num_blocks,
@@ -120,16 +121,20 @@ void SparseCholeskySolverUnitTest(
const int max_block_size,
const double block_density,
std::mt19937& prng) {
LinearSolver::Options sparse_cholesky_options;
sparse_cholesky_options.sparse_linear_algebra_library_type =
sparse_linear_algebra_library_type;
sparse_cholesky_options.ordering_type = ordering_type;
#ifndef CERES_NO_CUDSS
ContextImpl context;
sparse_cholesky_options.context = &context;
std::string error;
CHECK(context.InitCuda(&error)) << error;
#endif // CERES_NO_CUDSS
LinearSolver::Options sparse_cholesky_options;
sparse_cholesky_options.sparse_linear_algebra_library_type =
sparse_linear_algebra_library_type;
sparse_cholesky_options.ordering_type = ordering_type;
sparse_cholesky_options.max_num_refinement_iterations = 0;
sparse_cholesky_options.use_mixed_precision_solves = use_single_precision;
auto sparse_cholesky = SparseCholesky::Create(sparse_cholesky_options);
const CompressedRowSparseMatrix::StorageType storage_type =
sparse_cholesky->StorageType();
@@ -156,22 +161,30 @@ void SparseCholeskySolverUnitTest(
LinearSolverTerminationType::SUCCESS);
Matrix eigen_lhs;
lhs->ToDenseMatrix(&eigen_lhs);
EXPECT_NEAR((actual - expected).norm() / actual.norm(),
0.0,
std::numeric_limits<double>::epsilon() * 20)
const double kTolerance =
(use_single_precision ? std::numeric_limits<float>::epsilon()
: std::numeric_limits<double>::epsilon()) *
20;
EXPECT_NEAR((actual - expected).norm() / actual.norm(), 0.0, kTolerance)
<< "\n"
<< eigen_lhs;
}
// SparseLinearAlgebraLibraryType
// FLOAT/DOUBLE
// OrderingType
// BlockStructure
using Param =
::testing::tuple<SparseLinearAlgebraLibraryType, OrderingType, bool>;
::testing::tuple<SparseLinearAlgebraLibraryType, bool, OrderingType, bool>;
std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
Param param = info.param;
std::stringstream ss;
ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
<< ::testing::get<1>(param) << "_"
<< (::testing::get<2>(param) ? "UseBlockStructure" : "NoBlockStructure");
<< (::testing::get<1>(param) ? "FLOAT" : "DOUBLE") << "_"
<< ::testing::get<2>(param) << "_"
<< (::testing::get<3>(param) ? "UseBlockStructure" : "NoBlockStructure");
return ss.str();
}
@@ -198,6 +211,7 @@ TEST_P(SparseCholeskyTest, FactorAndSolve) {
SparseCholeskySolverUnitTest(::testing::get<0>(param),
::testing::get<1>(param),
::testing::get<2>(param),
::testing::get<3>(param),
num_blocks,
kMinBlockSize,
kMaxBlockSize,
@@ -214,38 +228,29 @@ INSTANTIATE_TEST_SUITE_P(
SuiteSparseCholesky,
SparseCholeskyTest,
::testing::Combine(::testing::Values(SUITE_SPARSE),
#if defined(CERES_NO_CHOLMOD_FLOAT)
::testing::Values(false),
#else
::testing::Values(false, true),
#endif // defined(CERES_NO_CHOLMOD_FLOAT)
#if defined(CERES_NO_CHOLMOD_PARTITION)
::testing::Values(OrderingType::AMD,
OrderingType::NATURAL),
#else
::testing::Values(OrderingType::AMD,
OrderingType::NESDIS,
OrderingType::NATURAL),
#endif // defined(CERES_NO_CHOLMOD_PARTITION)
::testing::Values(true, false)),
ParamInfoToString);
#endif
#if !defined(CERES_NO_SUITESPARSE) && !defined(CERES_NO_CHOLMOD_PARTITION)
INSTANTIATE_TEST_SUITE_P(
SuiteSparseCholeskyMETIS,
SparseCholeskyTest,
::testing::Combine(::testing::Values(SUITE_SPARSE),
::testing::Values(OrderingType::NESDIS),
::testing::Values(true, false)),
ParamInfoToString);
#endif // !defined(CERES_NO_SUITESPARSE) &&
// !defined(CERES_NO_CHOLMOD_PARTITION)
#endif // !defined(CERES_NO_SUITESPARSE)
#ifndef CERES_NO_ACCELERATE_SPARSE
INSTANTIATE_TEST_SUITE_P(
AccelerateSparseCholesky,
SparseCholeskyTest,
::testing::Combine(::testing::Values(ACCELERATE_SPARSE),
::testing::Values(OrderingType::AMD,
OrderingType::NESDIS,
OrderingType::NATURAL),
::testing::Values(true, false)),
ParamInfoToString);
INSTANTIATE_TEST_SUITE_P(
AccelerateSparseCholeskySingle,
SparseCholeskyTest,
::testing::Combine(::testing::Values(ACCELERATE_SPARSE),
::testing::Values(false, true),
::testing::Values(OrderingType::AMD,
OrderingType::NESDIS,
OrderingType::NATURAL),
@@ -258,59 +263,30 @@ INSTANTIATE_TEST_SUITE_P(
EigenSparseCholesky,
SparseCholeskyTest,
::testing::Combine(::testing::Values(EIGEN_SPARSE),
::testing::Values(false, true),
#if defined(CERES_NO_EIGEN_METIS)
::testing::Values(OrderingType::AMD,
OrderingType::NATURAL),
::testing::Values(true, false)),
ParamInfoToString);
INSTANTIATE_TEST_SUITE_P(
FloatEigenSparseCholesky,
SparseCholeskyTest,
::testing::Combine(::testing::Values(EIGEN_SPARSE),
#else
::testing::Values(OrderingType::AMD,
OrderingType::NATURAL),
OrderingType::NATURAL,
OrderingType::NESDIS),
#endif // defined(CERES_NO_EIGEN_METIS)
::testing::Values(true, false)),
ParamInfoToString);
#endif // CERES_USE_EIGEN_SPARSE
#ifndef CERES_NO_CUDSS
using CuDSSCholeskyDouble = CudaSparseCholesky<double>;
INSTANTIATE_TEST_SUITE_P(
CuDSSCholeskyDouble,
SparseCholeskyTest,
::testing::Combine(::testing::Values(CUDA_SPARSE),
::testing::Values(OrderingType::AMD),
::testing::Values(true, false)),
ParamInfoToString);
using CuDSSCholeskySingle = CudaSparseCholesky<float>;
INSTANTIATE_TEST_SUITE_P(
CuDSSCholeskySingle,
CudaCholesky,
SparseCholeskyTest,
::testing::Combine(::testing::Values(CUDA_SPARSE),
::testing::Values(false, true),
::testing::Values(OrderingType::AMD),
::testing::Values(true, false)),
ParamInfoToString);
#endif // CERES_NO_CUDSS
#if defined(CERES_USE_EIGEN_SPARSE) && !defined(CERES_NO_EIGEN_METIS)
INSTANTIATE_TEST_SUITE_P(
EigenSparseCholeskyMETIS,
SparseCholeskyTest,
::testing::Combine(::testing::Values(EIGEN_SPARSE),
::testing::Values(OrderingType::NESDIS),
::testing::Values(true, false)),
ParamInfoToString);
INSTANTIATE_TEST_SUITE_P(
EigenSparseCholeskySingleMETIS,
SparseCholeskyTest,
::testing::Combine(::testing::Values(EIGEN_SPARSE),
::testing::Values(OrderingType::NESDIS),
::testing::Values(true, false)),
ParamInfoToString);
#endif // defined(CERES_USE_EIGEN_SPARSE) && !defined(CERES_NO_EIGEN_METIS)
class MockSparseCholesky : public SparseCholesky {
public:
MOCK_CONST_METHOD0(StorageType, CompressedRowSparseMatrix::StorageType());
internal/ceres/suitesparse.cc
+150 -2
View File
@@ -250,9 +250,9 @@ bool SuiteSparse::BlockOrdering(const cholmod_sparse* A,
block_matrix.i = reinterpret_cast<void*>(block_rows.data());
block_matrix.x = nullptr;
block_matrix.stype = A->stype;
block_matrix.itype = CHOLMOD_INT;
block_matrix.itype = A->itype;
block_matrix.xtype = CHOLMOD_PATTERN;
block_matrix.dtype = CHOLMOD_DOUBLE;
block_matrix.dtype = A->dtype;
block_matrix.sorted = 1;
block_matrix.packed = 1;
@@ -466,6 +466,154 @@ LinearSolverTerminationType SuiteSparseCholesky::Solve(const double* rhs,
return LinearSolverTerminationType::SUCCESS;
}
#ifndef CERES_NO_CHOLMOD_FLOAT
std::unique_ptr<SparseCholesky> FloatSuiteSparseCholesky::Create(
const OrderingType ordering_type) {
return std::unique_ptr<SparseCholesky>(
new FloatSuiteSparseCholesky(ordering_type));
}
FloatSuiteSparseCholesky::FloatSuiteSparseCholesky(
const OrderingType ordering_type)
: ordering_type_(ordering_type), factor_(nullptr) {}
FloatSuiteSparseCholesky::~FloatSuiteSparseCholesky() {
if (factor_ != nullptr) {
ss_.Free(factor_);
}
}
LinearSolverTerminationType FloatSuiteSparseCholesky::Factorize(
CompressedRowSparseMatrix* lhs, std::string* message) {
if (lhs == nullptr) {
*message = "Failure: Input lhs is nullptr.";
return LinearSolverTerminationType::FATAL_ERROR;
}
cholmod_sparse cholmod_lhs = ss_.CreateSparseMatrixTransposeView(lhs);
float_lhs_values_ =
ConstVectorRef(lhs->values(), lhs->num_nonzeros()).cast<float>();
cholmod_lhs.dtype = CHOLMOD_SINGLE;
cholmod_lhs.x = reinterpret_cast<void*>(float_lhs_values_.data());
// 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 ||
(lhs->col_blocks().empty() || lhs->row_blocks().empty())) {
factor_ = ss_.AnalyzeCholesky(&cholmod_lhs, ordering_type_, message);
} else {
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);
}
CompressedRowSparseMatrix::StorageType FloatSuiteSparseCholesky::StorageType()
const {
return ((ordering_type_ == OrderingType::NATURAL)
? CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR
: CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
}
LinearSolverTerminationType FloatSuiteSparseCholesky::Solve(
const double* rhs, double* solution, std::string* message) {
// Error checking
if (factor_ == nullptr) {
*message = "Solve called without a call to Factorize first.";
return LinearSolverTerminationType::FATAL_ERROR;
}
const int num_cols = factor_->n;
cholmod_dense cholmod_rhs = ss_.CreateDenseVectorView(rhs, num_cols);
float_rhs_ = ConstVectorRef(rhs, num_cols).cast<float>();
cholmod_rhs.dtype = CHOLMOD_SINGLE;
cholmod_rhs.x = reinterpret_cast<void*>(float_rhs_.data());
cholmod_dense* cholmod_dense_solution =
ss_.Solve(factor_, &cholmod_rhs, message);
if (cholmod_dense_solution == nullptr) {
return LinearSolverTerminationType::FAILURE;
}
CHECK_EQ(cholmod_dense_solution->dtype, CHOLMOD_SINGLE);
VectorRef(solution, num_cols) =
Eigen::Map<Eigen::VectorXf>(
reinterpret_cast<float*>(cholmod_dense_solution->x), num_cols)
.cast<double>();
ss_.Free(cholmod_dense_solution);
return LinearSolverTerminationType::SUCCESS;
}
#else
std::unique_ptr<SparseCholesky> FloatSuiteSparseCholesky::Create(
const OrderingType ordering_type) {
return nullptr;
}
FloatSuiteSparseCholesky::FloatSuiteSparseCholesky(
const OrderingType ordering_type)
: ordering_type_(ordering_type), factor_(nullptr) {}
FloatSuiteSparseCholesky::~FloatSuiteSparseCholesky() {
if (factor_ != nullptr) {
ss_.Free(factor_);
}
}
LinearSolverTerminationType FloatSuiteSparseCholesky::Factorize(
CompressedRowSparseMatrix* lhs, std::string* message) {
*message =
"Single precision Cholesky factorization is not supported. If "
"you are seeing this failure, then you have discovered a Ceres "
"Solver bug. Please get in touch with the Ceres Solver team at "
"ceres-solver@googlegroups.com";
return LinearSolverTerminationType::FATAL_ERROR;
}
CompressedRowSparseMatrix::StorageType FloatSuiteSparseCholesky::StorageType()
const {
return ((ordering_type_ == OrderingType::NATURAL)
? CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR
: CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
}
LinearSolverTerminationType FloatSuiteSparseCholesky::Solve(
const double* rhs, double* solution, std::string* message) {
*message =
"Single precision Cholesky factorization is not supported. If "
"you are seeing this failure, then you have discovered a Ceres "
"Solver bug. Please get in touch with the Ceres Solver team at "
"ceres-solver@googlegroups.com";
return LinearSolverTerminationType::FATAL_ERROR;
}
#endif // CERES_NO_CHOLMOD_FLOAT
} // namespace ceres::internal
#endif // CERES_NO_SUITESPARSE
internal/ceres/suitesparse.h
+24 -1
View File
@@ -264,7 +264,30 @@ class CERES_NO_EXPORT SuiteSparseCholesky final : public SparseCholesky {
const OrderingType ordering_type_;
SuiteSparse ss_;
cholmod_factor* factor_;
cholmod_factor* factor_ = nullptr;
};
class CERES_NO_EXPORT FloatSuiteSparseCholesky final : public SparseCholesky {
public:
static std::unique_ptr<SparseCholesky> Create(OrderingType ordering_type);
// SparseCholesky interface.
~FloatSuiteSparseCholesky() override;
CompressedRowSparseMatrix::StorageType StorageType() const final;
LinearSolverTerminationType Factorize(CompressedRowSparseMatrix* lhs,
std::string* message) final;
LinearSolverTerminationType Solve(const double* rhs,
double* solution,
std::string* message) final;
private:
explicit FloatSuiteSparseCholesky(const OrderingType ordering_type);
const OrderingType ordering_type_;
SuiteSparse ss_;
cholmod_factor* factor_ = nullptr;
Eigen::VectorXf float_lhs_values_;
Eigen::VectorXf float_rhs_;
};
} // namespace ceres::internal