diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 27ea80e40..b05c3a29d 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -81,8 +81,9 @@ DEFINE_string(blocks_for_inner_iterations, "automatic", "Options are: " "automatic, cameras, points, cameras,points, points,cameras"); DEFINE_string(linear_solver, "sparse_schur", "Options are: " - "sparse_schur, dense_schur, iterative_schur, sparse_normal_cholesky, " - "dense_qr, dense_normal_cholesky and cgnr."); + "sparse_schur, dense_schur, iterative_schur, " + "sparse_normal_cholesky, dense_qr, dense_normal_cholesky, " + "and cgnr."); DEFINE_bool(explicit_schur_complement, false, "If using ITERATIVE_SCHUR " "then explicitly compute the Schur complement."); DEFINE_string(preconditioner, "jacobi", "Options are: " @@ -94,11 +95,13 @@ DEFINE_bool(use_spse_initialization, false, "Use power series expansion to initialize the solution in ITERATIVE_SCHUR linear solver."); DEFINE_string(sparse_linear_algebra_library, "suite_sparse", - "Options are: suite_sparse, cx_sparse, accelerate_sparse and eigen_sparse."); + "Options are: suite_sparse, accelerate_sparse, eigen_sparse, and " + "cuda_sparse."); DEFINE_string(dense_linear_algebra_library, "eigen", "Options are: eigen, lapack, and cuda"); DEFINE_string(ordering_type, "amd", "Options are: amd, nesdis"); -DEFINE_string(linear_solver_ordering, "automatic", "Options are: automatic and user"); +DEFINE_string(linear_solver_ordering, "automatic", + "Options are: automatic and user"); DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent " "rotations. If false, angle axis is used."); @@ -351,6 +354,7 @@ void SolveProblem(const char* filename) { SetSolverOptionsFromFlags(&bal_problem, &options); options.gradient_tolerance = 1e-16; options.function_tolerance = 1e-16; + options.parameter_tolerance = 1e-16; Solver::Summary summary; Solve(options, &problem, &summary); std::cout << summary.FullReport() << "\n"; diff --git a/include/ceres/types.h b/include/ceres/types.h index 1b07ab37f..0a110e252 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -137,7 +137,7 @@ enum PreconditionerType { // well the matrix Q approximates J'J, or how well the chosen // residual blocks approximate the non-linear least squares // problem. - SUBSET, + SUBSET }; enum VisibilityClusteringType { @@ -175,6 +175,9 @@ enum SparseLinearAlgebraLibraryType { // Apple's Accelerate framework sparse linear algebra routines. ACCELERATE_SPARSE, + // Nvidia's cuSPARSE library. + CUDA_SPARSE, + // No sparse linear solver should be used. This does not necessarily // imply that Ceres was built without any sparse library, although that // is the likely use case, merely that one should not be used. diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc index 99d53755b..f4f9e0fc2 100644 --- a/internal/ceres/cgnr_solver.cc +++ b/internal/ceres/cgnr_solver.cc @@ -35,12 +35,15 @@ #include "ceres/block_jacobi_preconditioner.h" #include "ceres/conjugate_gradients_solver.h" +#include "ceres/cuda_sparse_matrix.h" +#include "ceres/cuda_vector.h" #include "ceres/internal/eigen.h" #include "ceres/linear_solver.h" #include "ceres/subset_preconditioner.h" #include "ceres/wall_time.h" #include "glog/logging.h" + namespace ceres::internal { // A linear operator which takes a matrix A and a diagonal vector D and @@ -117,7 +120,14 @@ CgnrSolver::CgnrSolver(LinearSolver::Options options) } } -CgnrSolver::~CgnrSolver() = default; +CgnrSolver::~CgnrSolver() { + for (int i = 0; i < 4; ++i) { + if (scratch_[i]) { + delete scratch_[i]; + scratch_[i] = nullptr; + } + } +} LinearSolver::Summary CgnrSolver::SolveImpl( BlockSparseMatrix* A, @@ -163,10 +173,13 @@ LinearSolver::Summary CgnrSolver::SolveImpl( cg_solution_ = Vector::Zero(A->num_cols()); for (int i = 0; i < 4; ++i) { - scratch_[i] = Vector::Zero(A->num_cols()); + if (scratch_[i] == nullptr) { + scratch_[i] = new Vector(A->num_cols()); + } } event_logger.AddEvent("Setup"); + LinearOperatorAdapter preconditioner(*preconditioner_); auto summary = ConjugateGradientsSolver( cg_options, lhs, rhs, preconditioner, scratch_, cg_solution_); @@ -175,4 +188,135 @@ LinearSolver::Summary CgnrSolver::SolveImpl( return summary; } +#ifndef CERES_NO_CUDA + +// A linear operator which takes a matrix A and a diagonal vector D and +// performs products of the form +// +// (A^T A + D^T D)x +// +// This is used to implement iterative general sparse linear solving with +// conjugate gradients, where A is the Jacobian and D is a regularizing +// parameter. A brief proof is included in cgnr_linear_operator.h. +class CERES_NO_EXPORT CudaCgnrLinearOperator final : public + ConjugateGradientsLinearOperator { + public: + CudaCgnrLinearOperator(CudaSparseMatrix& A, + const CudaVector& D, + CudaVector* z) : A_(A), D_(D), z_(z) {} + + void RightMultiplyAndAccumulate(const CudaVector& x, CudaVector& y) final { + // z = Ax + z_->SetZero(); + A_.RightMultiplyAndAccumulate(x, z_); + + // y = y + Atz + // = y + AtAx + A_.LeftMultiplyAndAccumulate(*z_, &y); + + // y = y + DtDx + y.DtDxpy(D_, x); + } + + private: + CudaSparseMatrix& A_; + const CudaVector& D_; + CudaVector* z_ = nullptr; +}; + +class CERES_NO_EXPORT CudaIdentityPreconditioner final : public + ConjugateGradientsLinearOperator { + public: + void RightMultiplyAndAccumulate(const CudaVector& x, CudaVector& y) final { + y.Axpby(1.0, x, 1.0); + } +}; + +CudaCgnrSolver::CudaCgnrSolver(LinearSolver::Options options) : + options_(std::move(options)) {} + +CudaCgnrSolver::~CudaCgnrSolver() { + for (int i = 0; i < 4; ++i) { + if (scratch_[i]) { + delete scratch_[i]; + scratch_[i] = nullptr; + } + } +} + +std::unique_ptr CudaCgnrSolver::Create( + LinearSolver::Options options, std::string* error) { + CHECK(error != nullptr); + if (options.preconditioner_type != IDENTITY) { + *error = "CudaCgnrSolver does not support preconditioner type " + + std::string(PreconditionerTypeToString(options.preconditioner_type)) + ". "; + return nullptr; + } + CHECK(options.context->IsCUDAInitialized()) + << "CudaCgnrSolver requires CUDA initialization."; + auto solver = std::make_unique(options); + return solver; +} + +void CudaCgnrSolver::CpuToGpuTransfer( + const CompressedRowSparseMatrix& A, const double* b, const double* D) { + if (A_ == nullptr) { + // Assume structure is not cached, do an initialization and structural copy. + A_ = std::make_unique(options_.context, A); + b_ = std::make_unique(options_.context, A.num_rows()); + x_ = std::make_unique(options_.context, A.num_cols()); + Atb_ = std::make_unique(options_.context, A.num_cols()); + Ax_ = std::make_unique(options_.context, A.num_rows()); + D_ = std::make_unique(options_.context, A.num_cols()); + for (int i = 0; i < 4; ++i) { + scratch_[i] = new CudaVector(options_.context, A.num_cols()); + } + } else { + // Assume structure is cached, do a value copy. + A_->CopyValuesFromCpu(A); + } + b_->CopyFromCpu(ConstVectorRef(b, A.num_rows())); + D_->CopyFromCpu(ConstVectorRef(D, A.num_cols())); +} + +LinearSolver::Summary CudaCgnrSolver::SolveImpl( + CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double* x) { + EventLogger event_logger("CudaCgnrSolver::Solve"); + LinearSolver::Summary summary; + summary.num_iterations = 0; + summary.termination_type = LinearSolverTerminationType::FATAL_ERROR; + + CpuToGpuTransfer(*A, b, per_solve_options.D); + event_logger.AddEvent("CPU to GPU Transfer"); + + // Form z = Atb. + Atb_->SetZero(); + A_->LeftMultiplyAndAccumulate(*b_, Atb_.get()); + + // Solve (AtA + DtD)x = z (= Atb). + x_->SetZero(); + CudaCgnrLinearOperator lhs(*A_, *D_, Ax_.get()); + + event_logger.AddEvent("Setup"); + + ConjugateGradientsSolverOptions cg_options; + cg_options.min_num_iterations = options_.min_num_iterations; + cg_options.max_num_iterations = options_.max_num_iterations; + cg_options.residual_reset_period = options_.residual_reset_period; + cg_options.q_tolerance = per_solve_options.q_tolerance; + cg_options.r_tolerance = per_solve_options.r_tolerance; + + CudaIdentityPreconditioner preconditioner; + summary = ConjugateGradientsSolver( + cg_options, lhs, *Atb_, preconditioner, scratch_, *x_); + x_->CopyTo(x); + event_logger.AddEvent("Solve"); + return summary; +} + +#endif // CERES_NO_CUDA + } // namespace ceres::internal diff --git a/internal/ceres/cgnr_solver.h b/internal/ceres/cgnr_solver.h index 698229658..a971c8180 100644 --- a/internal/ceres/cgnr_solver.h +++ b/internal/ceres/cgnr_solver.h @@ -33,6 +33,8 @@ #include +#include "ceres/cuda_sparse_matrix.h" +#include "ceres/cuda_vector.h" #include "ceres/internal/export.h" #include "ceres/linear_solver.h" @@ -65,9 +67,40 @@ class CERES_NO_EXPORT CgnrSolver final : public BlockSparseMatrixSolver { const LinearSolver::Options options_; std::unique_ptr preconditioner_; Vector cg_solution_; - Vector scratch_[4]; + Vector* scratch_[4] = {nullptr, nullptr, nullptr, nullptr}; }; +#ifndef CERES_NO_CUDA +// A Cuda-accelerated version of CgnrSolver. +// This solver assumes that the sparsity structure of A remains constant for its +// lifetime. +class CERES_NO_EXPORT CudaCgnrSolver final : public CompressedRowSparseMatrixSolver { + public: + explicit CudaCgnrSolver(LinearSolver::Options options); + static std::unique_ptr Create( + LinearSolver::Options options, std::string* error); + ~CudaCgnrSolver() override; + + Summary SolveImpl(CompressedRowSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double* x) final; + + private: + void CpuToGpuTransfer( + const CompressedRowSparseMatrix& A, const double* b, const double* D); + + LinearSolver::Options options_; + std::unique_ptr A_; + std::unique_ptr b_; + std::unique_ptr x_; + std::unique_ptr Atb_; + std::unique_ptr Ax_; + std::unique_ptr D_; + CudaVector* scratch_[4] = {nullptr, nullptr, nullptr, nullptr}; +}; +#endif // CERES_NO_CUDA + } // namespace ceres::internal #endif // CERES_INTERNAL_CGNR_SOLVER_H_ diff --git a/internal/ceres/conjugate_gradients_solver.h b/internal/ceres/conjugate_gradients_solver.h index f2eea8745..0c3581b9c 100644 --- a/internal/ceres/conjugate_gradients_solver.h +++ b/internal/ceres/conjugate_gradients_solver.h @@ -100,25 +100,25 @@ struct ConjugateGradientsSolverOptions { // allows us to have a single implementation that works on CPU and GPU based // matrices and vectors. // -// scratch must contain four DenseVector objects of the same size as rhs and -// solution. By asking the user for scratch space, we guarantee that we will not -// perform any allocations inside this function. +// scratch must contain pointers to four DenseVector objects of the same size as +// rhs and solution. By asking the user for scratch space, we guarantee that we +// will not perform any allocations inside this function. template LinearSolver::Summary ConjugateGradientsSolver( const ConjugateGradientsSolverOptions options, ConjugateGradientsLinearOperator& lhs, const DenseVectorType& rhs, ConjugateGradientsLinearOperator& preconditioner, - DenseVectorType scratch[4], + DenseVectorType* scratch[4], DenseVectorType& solution) { auto IsZeroOrInfinity = [](double x) { return ((x == 0.0) || std::isinf(x)); }; - DenseVectorType& p = scratch[0]; - DenseVectorType& r = scratch[1]; - DenseVectorType& z = scratch[2]; - DenseVectorType& tmp = scratch[3]; + DenseVectorType& p = *scratch[0]; + DenseVectorType& r = *scratch[1]; + DenseVectorType& z = *scratch[2]; + DenseVectorType& tmp = *scratch[3]; LinearSolver::Summary summary; summary.termination_type = LinearSolverTerminationType::NO_CONVERGENCE; @@ -160,7 +160,7 @@ LinearSolver::Summary ConjugateGradientsSolver( SetZero(z); preconditioner.RightMultiplyAndAccumulate(r, z); - double last_rho = rho; + const double last_rho = rho; // rho = r.dot(z); rho = Dot(r, z); if (IsZeroOrInfinity(rho)) { @@ -172,7 +172,7 @@ LinearSolver::Summary ConjugateGradientsSolver( if (summary.num_iterations == 1) { Copy(z, p); } else { - double beta = rho / last_rho; + const double beta = rho / last_rho; if (IsZeroOrInfinity(beta)) { summary.termination_type = LinearSolverTerminationType::FAILURE; summary.message = StringPrintf( diff --git a/internal/ceres/conjugate_gradients_solver_test.cc b/internal/ceres/conjugate_gradients_solver_test.cc index a01dfc948..e59eba70f 100644 --- a/internal/ceres/conjugate_gradients_solver_test.cc +++ b/internal/ceres/conjugate_gradients_solver_test.cc @@ -74,9 +74,10 @@ TEST(ConjugateGradientTest, Solves3x3IdentitySystem) { IdentityPreconditioner identity(A->num_cols()); LinearOperatorAdapter lhs(*A); LinearOperatorAdapter preconditioner(identity); - - auto summary = - ConjugateGradientsSolver(cg_options, lhs, b, preconditioner, scratch, x); + Vector* scratch_array[4] = + {&scratch[0], &scratch[1], &scratch[2], &scratch[3]}; + auto summary = ConjugateGradientsSolver( + cg_options, lhs, b, preconditioner, scratch_array, x); EXPECT_EQ(summary.termination_type, LinearSolverTerminationType::SUCCESS); ASSERT_EQ(summary.num_iterations, 1); @@ -135,13 +136,14 @@ TEST(ConjuateGradientTest, Solves3x3SymmetricSystem) { for (int i = 0; i < 4; ++i) { scratch[i] = Vector::Zero(A->num_cols()); } - + Vector* scratch_array[4] = + {&scratch[0], &scratch[1], &scratch[2], &scratch[3]}; IdentityPreconditioner identity(A->num_cols()); LinearOperatorAdapter lhs(*A); LinearOperatorAdapter preconditioner(identity); - auto summary = - ConjugateGradientsSolver(cg_options, lhs, b, preconditioner, scratch, x); + auto summary = ConjugateGradientsSolver( + cg_options, lhs, b, preconditioner, scratch_array, x); EXPECT_EQ(summary.termination_type, LinearSolverTerminationType::SUCCESS); diff --git a/internal/ceres/cuda_dense_cholesky_test.cc b/internal/ceres/cuda_dense_cholesky_test.cc index c7b11ce4f..148392313 100644 --- a/internal/ceres/cuda_dense_cholesky_test.cc +++ b/internal/ceres/cuda_dense_cholesky_test.cc @@ -36,8 +36,7 @@ #include "glog/logging.h" #include "gtest/gtest.h" -namespace ceres { -namespace internal { +namespace ceres::internal { #ifndef CERES_NO_CUDA @@ -45,6 +44,8 @@ TEST(CUDADenseCholesky, InvalidOptionOnCreate) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; auto dense_cuda_solver = CUDADenseCholesky::Create(options); EXPECT_EQ(dense_cuda_solver, nullptr); } @@ -63,6 +64,8 @@ TEST(CUDADenseCholesky, Cholesky4x4Matrix) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseCholesky::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -90,6 +93,8 @@ TEST(CUDADenseCholesky, SingularMatrix) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseCholesky::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -110,6 +115,8 @@ TEST(CUDADenseCholesky, NegativeMatrix) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseCholesky::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -123,6 +130,8 @@ TEST(CUDADenseCholesky, MustFactorizeBeforeSolve) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseCholesky::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -140,6 +149,8 @@ TEST(CUDADenseCholesky, Randomized1600x1600Tests) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = ceres::CUDA; std::unique_ptr dense_cholesky = CUDADenseCholesky::Create(options); @@ -175,12 +186,20 @@ TEST(CUDADenseCholeskyMixedPrecision, InvalidOptionsOnCreate) { { // Did not ask for CUDA, and did not ask for mixed precision. LinearSolver::Options options; + ContextImpl context; + options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; auto solver = CUDADenseCholeskyMixedPrecision::Create(options); ASSERT_EQ(solver, nullptr); } { // Asked for CUDA, but did not ask for mixed precision. LinearSolver::Options options; + ContextImpl context; + options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = ceres::CUDA; auto solver = CUDADenseCholeskyMixedPrecision::Create(options); ASSERT_EQ(solver, nullptr); @@ -204,6 +223,8 @@ TEST(CUDADenseCholeskyMixedPrecision, Cholesky4x4Matrix1Step) { options.max_num_refinement_iterations = 0; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; options.use_mixed_precision_solves = true; auto solver = CUDADenseCholeskyMixedPrecision::Create(options); @@ -239,6 +260,8 @@ TEST(CUDADenseCholeskyMixedPrecision, Cholesky4x4Matrix4Steps) { options.max_num_refinement_iterations = 3; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; options.use_mixed_precision_solves = true; auto solver = CUDADenseCholeskyMixedPrecision::Create(options); @@ -267,6 +290,8 @@ TEST(CUDADenseCholeskyMixedPrecision, Randomized1600x1600Tests) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = ceres::CUDA; options.use_mixed_precision_solves = true; options.max_num_refinement_iterations = 20; @@ -302,5 +327,4 @@ TEST(CUDADenseCholeskyMixedPrecision, Randomized1600x1600Tests) { #endif // CERES_NO_CUDA -} // namespace internal -} // namespace ceres +} // namespace ceres::internal diff --git a/internal/ceres/cuda_dense_qr_test.cc b/internal/ceres/cuda_dense_qr_test.cc index 9eb5d4d58..cc0dbd6d0 100644 --- a/internal/ceres/cuda_dense_qr_test.cc +++ b/internal/ceres/cuda_dense_qr_test.cc @@ -41,6 +41,10 @@ namespace ceres::internal { TEST(CUDADenseQR, InvalidOptionOnCreate) { LinearSolver::Options options; + ContextImpl context; + options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; auto dense_cuda_solver = CUDADenseQR::Create(options); EXPECT_EQ(dense_cuda_solver, nullptr); } @@ -58,6 +62,8 @@ TEST(CUDADenseQR, QR4x4Matrix) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseQR::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -90,6 +96,8 @@ TEST(CUDADenseQR, QR4x2Matrix) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseQR::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -112,6 +120,8 @@ TEST(CUDADenseQR, MustFactorizeBeforeSolve) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = CUDA; auto dense_cuda_solver = CUDADenseQR::Create(options); ASSERT_NE(dense_cuda_solver, nullptr); @@ -130,10 +140,12 @@ TEST(CUDADenseQR, Randomized1600x100Tests) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + EXPECT_TRUE(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = ceres::CUDA; std::unique_ptr dense_qr = CUDADenseQR::Create(options); - const int kNumTrials = 100; + const int kNumTrials = 20; for (int i = 0; i < kNumTrials; ++i) { LhsType lhs = LhsType::Random(kNumRows, kNumCols); SolutionType x_expected = SolutionType::Random(kNumCols); diff --git a/internal/ceres/dense_cholesky.cc b/internal/ceres/dense_cholesky.cc index 244de7426..af42e73bf 100644 --- a/internal/ceres/dense_cholesky.cc +++ b/internal/ceres/dense_cholesky.cc @@ -348,9 +348,8 @@ LinearSolverTerminationType RefinedDenseCholesky::Solve(const double* rhs, #ifndef CERES_NO_CUDA bool CUDADenseCholesky::Init(ContextImpl* context, std::string* message) { - if (!context->InitCUDA(message)) { - return false; - } + CHECK(context->IsCUDAInitialized()) + << "CUDADenseCholesky requires CUDA initialization."; cusolver_handle_ = context->cusolver_handle_; stream_ = context->stream_; error_.Reserve(1); @@ -491,9 +490,8 @@ CUDADenseCholeskyMixedPrecision::Create(const LinearSolver::Options& options) { bool CUDADenseCholeskyMixedPrecision::Init(const LinearSolver::Options& options, std::string* message) { - if (!options.context->InitCUDA(message)) { - return false; - } + CHECK(options.context->IsCUDAInitialized()) + << "CUDADenseCholeskyMixedPrecision requires CUDA initialization."; cusolver_handle_ = options.context->cusolver_handle_; cublas_handle_ = options.context->cublas_handle_; stream_ = options.context->stream_; diff --git a/internal/ceres/dense_cholesky_test.cc b/internal/ceres/dense_cholesky_test.cc index 1c0df568d..34ff4f38f 100644 --- a/internal/ceres/dense_cholesky_test.cc +++ b/internal/ceres/dense_cholesky_test.cc @@ -76,6 +76,8 @@ TEST_P(DenseCholeskyTest, FactorAndSolve) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + CHECK(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = ::testing::get<0>(GetParam()); options.use_mixed_precision_solves = ::testing::get<1>(GetParam()); const int kNumRefinementSteps = 4; diff --git a/internal/ceres/dense_qr_test.cc b/internal/ceres/dense_qr_test.cc index acef75070..f562d972b 100644 --- a/internal/ceres/dense_qr_test.cc +++ b/internal/ceres/dense_qr_test.cc @@ -68,6 +68,8 @@ TEST_P(DenseQRTest, FactorAndSolve) { LinearSolver::Options options; ContextImpl context; options.context = &context; + std::string error; + CHECK(context.InitCUDA(&error)) << error; options.dense_linear_algebra_library_type = GetParam(); const double kEpsilon = std::numeric_limits::epsilon() * 1.5e4; std::unique_ptr dense_qr = DenseQR::Create(options); diff --git a/internal/ceres/evaluator.cc b/internal/ceres/evaluator.cc index f44796a38..91cf22104 100644 --- a/internal/ceres/evaluator.cc +++ b/internal/ceres/evaluator.cc @@ -64,10 +64,17 @@ std::unique_ptr Evaluator::Create(const Evaluator::Options& options, case DENSE_SCHUR: case SPARSE_SCHUR: case ITERATIVE_SCHUR: - case CGNR: - return std::make_unique< - ProgramEvaluator>( - options, program); + case CGNR: { + if (options.sparse_linear_algebra_library_type == CUDA_SPARSE) { + return std::make_unique< + ProgramEvaluator>(options, program); + } else { + return std::make_unique< + ProgramEvaluator>( + options, program); + } + } case SPARSE_NORMAL_CHOLESKY: if (options.dynamic_sparsity) { return std::make_unique< diff --git a/internal/ceres/evaluator.h b/internal/ceres/evaluator.h index 68a4fb28a..7ea523c58 100644 --- a/internal/ceres/evaluator.h +++ b/internal/ceres/evaluator.h @@ -65,6 +65,7 @@ class CERES_NO_EXPORT Evaluator { int num_threads = 1; int num_eliminate_blocks = -1; LinearSolverType linear_solver_type = DENSE_QR; + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = NO_SPARSE; bool dynamic_sparsity = false; ContextImpl* context = nullptr; EvaluationCallback* evaluation_callback = nullptr; diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc index 607502e72..0463afbf3 100644 --- a/internal/ceres/iterative_schur_complement_solver.cc +++ b/internal/ceres/iterative_schur_complement_solver.cc @@ -129,6 +129,9 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl( for (int i = 0; i < 4; ++i) { scratch[i] = Vector::Zero(schur_complement_->num_cols()); } + Vector* scratch_ptr[4] = { + &scratch[0], &scratch[1], &scratch[2], &scratch[3] + }; event_logger.AddEvent("Setup"); @@ -137,7 +140,7 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl( lhs, schur_complement_->rhs(), preconditioner, - scratch, + scratch_ptr, reduced_linear_system_solution_); if (summary.termination_type != LinearSolverTerminationType::FAILURE && diff --git a/internal/ceres/linear_solver.cc b/internal/ceres/linear_solver.cc index d9ef15006..2e3fe4be6 100644 --- a/internal/ceres/linear_solver.cc +++ b/internal/ceres/linear_solver.cc @@ -76,8 +76,15 @@ std::unique_ptr LinearSolver::Create( CHECK(options.context != nullptr); switch (options.type) { - case CGNR: + case CGNR: { +#ifndef CERES_NO_CUDA + if (options.sparse_linear_algebra_library_type == CUDA_SPARSE) { + std::string error; + return CudaCgnrSolver::Create(options, &error); + } +#endif return std::make_unique(options); + } break; case SPARSE_NORMAL_CHOLESKY: #if defined(CERES_NO_SPARSE) diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index ba96e1cc5..631862ced 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc @@ -217,7 +217,14 @@ SparseSchurComplementSolver::SparseSchurComplementSolver( } } -SparseSchurComplementSolver::~SparseSchurComplementSolver() = default; +SparseSchurComplementSolver::~SparseSchurComplementSolver() { + for (int i = 0; i < 4; ++i) { + if (scratch_[i]) { + delete scratch_[i]; + scratch_[i] = nullptr; + } + } +} // Determine the non-zero blocks in the Schur Complement matrix, and // initialize a BlockRandomAccessSparseMatrix object. @@ -394,11 +401,11 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingConjugateGradients( cg_options.r_tolerance = per_solve_options.r_tolerance; cg_solution_ = Vector::Zero(sc->num_rows()); - Vector scratch[4]; for (int i = 0; i < 4; ++i) { - scratch_[i] = Vector::Zero(sc->num_rows()); + if (scratch_[i] == nullptr) { + scratch_[i] = new Vector(sc->num_rows()); + } } - auto summary = ConjugateGradientsSolver( cg_options, *lhs, rhs(), *preconditioner, scratch_, cg_solution_); VectorRef(solution, sc->num_rows()) = cg_solution_; diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h index 4fde69cb5..cac96b5ba 100644 --- a/internal/ceres/schur_complement_solver.h +++ b/internal/ceres/schur_complement_solver.h @@ -188,7 +188,7 @@ class CERES_NO_EXPORT SparseSchurComplementSolver final std::unique_ptr sparse_cholesky_; std::unique_ptr preconditioner_; Vector cg_solution_; - Vector scratch_[4]; + Vector* scratch_[4] = {nullptr, nullptr, nullptr, nullptr}; }; } // namespace ceres::internal diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index b772f3f15..0e0790e87 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc @@ -327,7 +327,8 @@ bool OptionsAreValidForCgnr(const Solver::Options& options, string* error) { return false; } - if (options.preconditioner_type == SUBSET) { + if (options.sparse_linear_algebra_library_type != CUDA_SPARSE && + options.preconditioner_type == SUBSET) { if (options.residual_blocks_for_subset_preconditioner.empty()) { *error = "When using SUBSET preconditioner, " @@ -341,6 +342,20 @@ bool OptionsAreValidForCgnr(const Solver::Options& options, string* error) { } } + // Check options for CGNR with CUDA_SPARSE. + if (options.sparse_linear_algebra_library_type == CUDA_SPARSE) { + if (!IsSparseLinearAlgebraLibraryTypeAvailable(CUDA_SPARSE)) { + *error = "Can't use CGNR with sparse_linear_algebra_library_type = " + "CUDA_SPARSE because support was not enabled when Ceres was built."; + return false; + } + if (options.preconditioner_type != IDENTITY) { + StringPrintf("Can't use CGNR with preconditioner_type = %s when " + "sparse_linear_algebra_library_type = CUDA_SPARSE.", + PreconditionerTypeToString(options.preconditioner_type)); + return false; + } + } return true; } @@ -652,6 +667,18 @@ std::string SchurStructureToString(const int row_block_size, return internal::StringPrintf("%s,%s,%s", row.c_str(), e.c_str(), f.c_str()); } +bool IsCudaRequired(const Solver::Options& options) { + if (options.linear_solver_type == DENSE_NORMAL_CHOLESKY || + options.linear_solver_type == DENSE_SCHUR || + options.linear_solver_type == DENSE_QR) { + return (options.dense_linear_algebra_library_type == CUDA); + } + if (options.linear_solver_type == CGNR) { + return (options.sparse_linear_algebra_library_type == CUDA_SPARSE); + } + return false; +} + } // namespace bool Solver::Options::IsValid(string* error) const { @@ -697,6 +724,13 @@ void Solver::Solve(const Solver::Options& options, Program* program = problem_impl->mutable_program(); PreSolveSummarize(options, problem_impl, summary); + if (IsCudaRequired(options)) { + if (!problem_impl->context()->InitCUDA(&summary->message)) { + LOG(ERROR) << "Terminating: " << summary->message; + return; + } + } + // If gradient_checking is enabled, wrap all cost functions in a // gradient checker and install a callback that terminates if any gradient // error is detected. @@ -866,22 +900,40 @@ string Solver::Summary::FullReport() const { } } - if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY || + const bool used_sparse_linear_algebra_library = + linear_solver_type_used == SPARSE_NORMAL_CHOLESKY || linear_solver_type_used == SPARSE_SCHUR || - (linear_solver_type_used == CGNR && - preconditioner_type_used == SUBSET) || + linear_solver_type_used == CGNR || (linear_solver_type_used == ITERATIVE_SCHUR && - (preconditioner_type_used == CLUSTER_JACOBI || - preconditioner_type_used == CLUSTER_TRIDIAGONAL))) { + (preconditioner_type_used == CLUSTER_JACOBI || + preconditioner_type_used == CLUSTER_TRIDIAGONAL)); + + const bool linear_solver_ordering_required = + linear_solver_type_used == SPARSE_SCHUR || + (linear_solver_type_used == ITERATIVE_SCHUR && + (preconditioner_type_used == CLUSTER_JACOBI || + preconditioner_type_used == CLUSTER_TRIDIAGONAL)) || + (linear_solver_type_used == CGNR && preconditioner_type_used == SUBSET); + + if (used_sparse_linear_algebra_library) { const char* mixed_precision_suffix = (mixed_precision_solves_used ? "(Mixed Precision)" : ""); - StringAppendF( - &report, - "\nSparse linear algebra library %15s + %s %s\n", - SparseLinearAlgebraLibraryTypeToString( - sparse_linear_algebra_library_type), - LinearSolverOrderingTypeToString(linear_solver_ordering_type), - mixed_precision_suffix); + if (linear_solver_ordering_required) { + StringAppendF( + &report, + "\nSparse linear algebra library %15s + %s %s\n", + SparseLinearAlgebraLibraryTypeToString( + sparse_linear_algebra_library_type), + LinearSolverOrderingTypeToString(linear_solver_ordering_type), + mixed_precision_suffix); + } else { + StringAppendF( + &report, + "\nSparse linear algebra library %15s %s\n", + SparseLinearAlgebraLibraryTypeToString( + sparse_linear_algebra_library_type), + mixed_precision_suffix); + } } StringAppendF(&report, "\n"); diff --git a/internal/ceres/solver_test.cc b/internal/ceres/solver_test.cc index b4b2d3446..92c1ddc25 100644 --- a/internal/ceres/solver_test.cc +++ b/internal/ceres/solver_test.cc @@ -878,6 +878,20 @@ TEST(Solver, CgnrOptionsIdentityPreconditioner) { options.dynamic_sparsity = false; options.use_mixed_precision_solves = true; EXPECT_FALSE(options.IsValid(&message)); + + options.sparse_linear_algebra_library_type = CUDA_SPARSE; + options.dynamic_sparsity = false; + options.use_mixed_precision_solves = false; + EXPECT_EQ(options.IsValid(&message), + IsSparseLinearAlgebraLibraryTypeAvailable(CUDA_SPARSE)); + + options.dynamic_sparsity = true; + options.use_mixed_precision_solves = false; + EXPECT_FALSE(options.IsValid(&message)); + + options.dynamic_sparsity = false; + options.use_mixed_precision_solves = true; + EXPECT_FALSE(options.IsValid(&message)); } TEST(Solver, CgnrOptionsJacobiPreconditioner) { @@ -940,6 +954,20 @@ TEST(Solver, CgnrOptionsJacobiPreconditioner) { options.dynamic_sparsity = false; options.use_mixed_precision_solves = true; EXPECT_FALSE(options.IsValid(&message)); + + options.sparse_linear_algebra_library_type = CUDA_SPARSE; + + options.dynamic_sparsity = false; + options.use_mixed_precision_solves = false; + EXPECT_FALSE(options.IsValid(&message)); + + options.dynamic_sparsity = true; + options.use_mixed_precision_solves = false; + EXPECT_FALSE(options.IsValid(&message)); + + options.dynamic_sparsity = false; + options.use_mixed_precision_solves = true; + EXPECT_FALSE(options.IsValid(&message)); } TEST(Solver, CgnrOptionsSubsetPreconditioner) { @@ -1012,6 +1040,19 @@ TEST(Solver, CgnrOptionsSubsetPreconditioner) { options.use_mixed_precision_solves = true; EXPECT_FALSE(options.IsValid(&message)); } + + options.sparse_linear_algebra_library_type = CUDA_SPARSE; + options.dynamic_sparsity = false; + options.use_mixed_precision_solves = false; + EXPECT_FALSE(options.IsValid(&message)); + + options.dynamic_sparsity = true; + options.use_mixed_precision_solves = false; + EXPECT_FALSE(options.IsValid(&message)); + + options.dynamic_sparsity = false; + options.use_mixed_precision_solves = true; + EXPECT_FALSE(options.IsValid(&message)); } TEST(Solver, CgnrOptionsSchurPreconditioners) { diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc index d27a47201..84e783f5b 100644 --- a/internal/ceres/trust_region_preprocessor.cc +++ b/internal/ceres/trust_region_preprocessor.cc @@ -261,6 +261,8 @@ bool SetupEvaluator(PreprocessedProblem* pp) { const Solver::Options& options = pp->options; pp->evaluator_options = Evaluator::Options(); pp->evaluator_options.linear_solver_type = options.linear_solver_type; + pp->evaluator_options.sparse_linear_algebra_library_type = + options.sparse_linear_algebra_library_type; pp->evaluator_options.num_eliminate_blocks = 0; if (IsSchurType(options.linear_solver_type)) { pp->evaluator_options.num_eliminate_blocks = diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 1d514de77..3de91ed58 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -108,6 +108,7 @@ const char* SparseLinearAlgebraLibraryTypeToString( CASESTR(SUITE_SPARSE); CASESTR(EIGEN_SPARSE); CASESTR(ACCELERATE_SPARSE); + CASESTR(CUDA_SPARSE); CASESTR(NO_SPARSE); default: return "UNKNOWN"; @@ -120,6 +121,7 @@ bool StringToSparseLinearAlgebraLibraryType( STRENUM(SUITE_SPARSE); STRENUM(EIGEN_SPARSE); STRENUM(ACCELERATE_SPARSE); + STRENUM(CUDA_SPARSE); STRENUM(NO_SPARSE); return false; } @@ -420,6 +422,14 @@ bool IsSparseLinearAlgebraLibraryTypeAvailable( #endif } + if (type == CUDA_SPARSE) { +#ifdef CERES_NO_CUDA + return false; +#else + return true; +#endif + } + if (type == NO_SPARSE) { return true; }