diff --git a/CMakeLists.txt b/CMakeLists.txt index 529c99b77..f266e9ca7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/cmake/config.h.in b/cmake/config.h.in index 51db4a034..94120aeb2 100644 --- a/cmake/config.h.in +++ b/cmake/config.h.in @@ -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@ diff --git a/docs/source/nnls_solving.rst b/docs/source/nnls_solving.rst index 184b94b72..bc0e16a57 100644 --- a/docs/source/nnls_solving.rst +++ b/docs/source/nnls_solving.rst @@ -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: diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 68438a10b..b31880fc1 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -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 diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 533184309..ed30c7d56 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -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 diff --git a/internal/ceres/float_suitesparse.cc b/internal/ceres/float_suitesparse.cc deleted file mode 100644 index cea5c4acb..000000000 --- a/internal/ceres/float_suitesparse.cc +++ /dev/null @@ -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 - -#include "absl/log/log.h" - -#if !defined(CERES_NO_SUITESPARSE) - -namespace ceres::internal { - -std::unique_ptr FloatSuiteSparseCholesky::Create( - OrderingType ordering_type) { - LOG(FATAL) << "FloatSuiteSparseCholesky is not available."; - return {}; -} - -} // namespace ceres::internal - -#endif // !defined(CERES_NO_SUITESPARSE) diff --git a/internal/ceres/float_suitesparse.h b/internal/ceres/float_suitesparse.h deleted file mode 100644 index b9d298ea9..000000000 --- a/internal/ceres/float_suitesparse.h +++ /dev/null @@ -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 - -#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 Create(OrderingType ordering_type); -}; - -} // namespace ceres::internal - -#endif // !defined(CERES_NO_SUITESPARSE) - -#endif // CERES_INTERNAL_FLOAT_SUITESPARSE_H_ diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index bb82e6deb..010ad2602 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc @@ -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) { diff --git a/internal/ceres/solver_test.cc b/internal/ceres/solver_test.cc index 0a83a347b..f6fd9c7ac 100644 --- a/internal/ceres/solver_test.cc +++ b/internal/ceres/solver_test.cc @@ -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; diff --git a/internal/ceres/sparse_cholesky.cc b/internal/ceres/sparse_cholesky.cc index f1b2c8854..23396a9f6 100644 --- a/internal/ceres/sparse_cholesky.cc +++ b/internal/ceres/sparse_cholesky.cc @@ -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" diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc index 2b8e852af..b2499ff8c 100644 --- a/internal/ceres/sparse_cholesky_test.cc +++ b/internal/ceres/sparse_cholesky_test.cc @@ -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::epsilon() * 20) + const double kTolerance = + (use_single_precision ? std::numeric_limits::epsilon() + : std::numeric_limits::epsilon()) * + 20; + + EXPECT_NEAR((actual - expected).norm() / actual.norm(), 0.0, kTolerance) << "\n" << eigen_lhs; } +// SparseLinearAlgebraLibraryType +// FLOAT/DOUBLE +// OrderingType +// BlockStructure using Param = - ::testing::tuple; + ::testing::tuple; std::string ParamInfoToString(testing::TestParamInfo 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; INSTANTIATE_TEST_SUITE_P( - CuDSSCholeskyDouble, - SparseCholeskyTest, - ::testing::Combine(::testing::Values(CUDA_SPARSE), - ::testing::Values(OrderingType::AMD), - ::testing::Values(true, false)), - ParamInfoToString); - -using CuDSSCholeskySingle = CudaSparseCholesky; -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()); diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc index eeb850ccf..2cc9ff38f 100644 --- a/internal/ceres/suitesparse.cc +++ b/internal/ceres/suitesparse.cc @@ -250,9 +250,9 @@ bool SuiteSparse::BlockOrdering(const cholmod_sparse* A, block_matrix.i = reinterpret_cast(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 FloatSuiteSparseCholesky::Create( + const OrderingType ordering_type) { + return std::unique_ptr( + 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(); + cholmod_lhs.dtype = CHOLMOD_SINGLE; + cholmod_lhs.x = reinterpret_cast(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(); + cholmod_rhs.dtype = CHOLMOD_SINGLE; + cholmod_rhs.x = reinterpret_cast(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( + reinterpret_cast(cholmod_dense_solution->x), num_cols) + .cast(); + ss_.Free(cholmod_dense_solution); + return LinearSolverTerminationType::SUCCESS; +} + +#else + +std::unique_ptr 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 diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h index 57d1d78d8..11b136dbd 100644 --- a/internal/ceres/suitesparse.h +++ b/internal/ceres/suitesparse.h @@ -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 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