Add mixed precision solves for SUITE_SPARSE
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
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 529c99b..f266e9c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -374,6 +374,11 @@
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 51db4a0..94120ae 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 184b94b..bc0e16a 100644
--- a/docs/source/nnls_solving.rst
+++ b/docs/source/nnls_solving.rst
@@ -905,20 +905,12 @@
: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 68438a1..b31880f 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -569,23 +569,24 @@
// 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 5331843..ed30c7d 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -184,7 +184,6 @@
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 cea5c4a..0000000
--- 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 <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)
diff --git a/internal/ceres/float_suitesparse.h b/internal/ceres/float_suitesparse.h
deleted file mode 100644
index b9d298e..0000000
--- 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 <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_
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index bb82e6d..010ad26 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -195,11 +195,16 @@
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 0a83a34..f6fd9c7 100644
--- a/internal/ceres/solver_test.cc
+++ b/internal/ceres/solver_test.cc
@@ -511,7 +511,11 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 f1b2c88..23396a9 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 2b8e852..b2499ff 100644
--- a/internal/ceres/sparse_cholesky_test.cc
+++ b/internal/ceres/sparse_cholesky_test.cc
@@ -113,6 +113,7 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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,
+ CudaCholesky,
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,
- 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 eeb850c..2cc9ff3 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -250,9 +250,9 @@
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 @@
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
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 57d1d78..11b136d 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -264,7 +264,30 @@
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
