Add support for dense CUDA solvers #1

1. Add CUDADenseCholesky64Bit, CUDADenseCholesky32Bit, & tests.
   CUDADenseCholesky32Bit uses the legacy versions of potrf/potrs
   in cuSolverDN, while CUDADenseCholesky64Bit uses the new 64-bit
   versions available since Cuda 11.1. The legacy versions are
   provided since some platforms such as the Nvidia Jetsons only
   support Cuda 10.2.
2. Expose CUDA as a new option under DenseLinearAlgebraLibraryType.
   The relevant option to string and string to option helper functions
   are modified accordingly.
3. Add cuda as a dense_linear_algebra_library option in bundle_adjuster
   to demonstrate the use of the new CUDA option.

Change-Id: I23615e1d301df5185ed646b3e33ee802508dae86

CMakeLists.txt

@@ -154,6 +154,7 @@ if (APPLE)
   option(ACCELERATESPARSE
     "Enable use of sparse solvers in Apple's Accelerate framework." ON)
 endif()
+option(CUDA "Enable use of CUDA linear algebra solvers." ON)
 option(LAPACK "Enable use of LAPACK directly within Ceres." ON)
 # Template specializations for the Schur complement based solvers. If
 # compile time, binary size or compiler performance is an issue, you
@@ -260,6 +261,31 @@ if (Eigen3_FOUND)
   endif (EIGENSPARSE)
 endif (Eigen3_FOUND)
 
+if (CUDA)
+  find_package(CUDA QUIET)
+  if (CUDA_FOUND)
+    message("-- Found CUDA version ${CUDA_VERSION}: "
+        "${CUDA_LIBRARIES};"
+        "${CUDA_cusolver_LIBRARY};"
+        "${CUDA_cusparse_LIBRARY}")
+    if (CUDA_VERSION VERSION_LESS 11.1)
+      message("-- CUDA version ${CUDA_VERSION} is older than 11.1, "
+          "using 32-bit API.")
+      add_definitions(-DCERES_CUDA_NO_64BIT_SOLVER_API)
+    else()
+      message("-- CUDA version ${CUDA_VERSION} is at least 11.1, "
+          "using 64-bit API.")
+    endif()
+  else (CUDA_FOUND)
+    message("-- Did not find CUDA library, disabling CUDA support.")
+    update_cache_variable(CUDA OFF)
+    list(APPEND CERES_COMPILE_OPTIONS CERES_NO_CUDA)
+  endif (CUDA_FOUND)
+else (CUDA)
+  message("-- Building without CUDA.")
+  list(APPEND CERES_COMPILE_OPTIONS CERES_NO_CUDA)
+endif (CUDA)
+
 if (LAPACK)
   find_package(LAPACK QUIET)
   if (LAPACK_FOUND)

cmake/config.h.in

@@ -53,6 +53,9 @@
 // If defined, Ceres was compiled without CXSparse.
 @CERES_NO_CXSPARSE@
 
+// If defined, Ceres was compiled without CUDA.
+@CERES_NO_CUDA@
+
 // If defined, Ceres was compiled without Apple's Accelerate framework solvers.
 @CERES_NO_ACCELERATE_SPARSE@
 

examples/bundle_adjuster.cc

@@ -93,7 +93,7 @@ DEFINE_string(visibility_clustering, "canonical_views",
 DEFINE_string(sparse_linear_algebra_library, "suite_sparse",
               "Options are: suite_sparse and cx_sparse.");
 DEFINE_string(dense_linear_algebra_library, "eigen",
-              "Options are: eigen and lapack.");
+              "Options are: eigen, lapack, and cuda");
 DEFINE_string(ordering, "automatic", "Options are: automatic, user.");
 
 DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "

include/ceres/types.h

@@ -186,6 +186,7 @@ enum SparseLinearAlgebraLibraryType {
 enum DenseLinearAlgebraLibraryType {
   EIGEN,
   LAPACK,
+  CUDA,
 };
 
 // Logging options

internal/ceres/CMakeLists.txt

@@ -218,6 +218,14 @@ if (ACCELERATESPARSE AND AccelerateSparse_FOUND)
   list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${AccelerateSparse_LIBRARIES})
 endif()
 
+if (CUDA_FOUND)
+  list(APPEND
+       CERES_LIBRARY_PRIVATE_DEPENDENCIES
+       ${CUDA_LIBRARIES}
+       ${CUDA_cusolver_LIBRARY}
+       ${CUDA_cusparse_LIBRARY})
+endif (CUDA_FOUND)
+
 if (LAPACK_FOUND)
   list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${LAPACK_LIBRARIES})
 endif ()
@@ -362,6 +370,10 @@ if (ACCELERATESPARSE)
   list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES_INCLUDE_DIRS
     ${AccelerateSparse_INCLUDE_DIRS})
 endif()
+if (CUDA)
+  list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES_INCLUDE_DIRS
+    ${CUDA_INCLUDE_DIRS})
+endif()
 # Add include locations for optional dependencies to the Ceres target without
 # duplication.
 list(REMOVE_DUPLICATES CERES_LIBRARY_PRIVATE_DEPENDENCIES_INCLUDE_DIRS)
@@ -453,6 +465,7 @@ if (BUILD_TESTING AND GFLAGS)
   ceres_test(cost_function_to_functor)
   ceres_test(covariance)
   ceres_test(cubic_interpolation)
+  ceres_test(cuda_dense_cholesky)
   ceres_test(dense_linear_solver)
   ceres_test(dense_cholesky)
   ceres_test(dense_qr)

internal/ceres/cuda_buffer.h

@@ -0,0 +1,98 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 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: joydeepb@cs.utexas.edu (Joydeep Biswas)
+
+#ifndef CERES_INTERNAL_CUDA_BUFFER_H_
+#define CERES_INTERNAL_CUDA_BUFFER_H_
+
+#include "ceres/internal/config.h"
+
+#ifndef CERES_NO_CUDA
+
+#include <vector>
+
+#include "cuda_runtime.h"
+#include "glog/logging.h"
+
+// An encapsulated buffer to maintain GPU memory, and handle transfers between
+// GPU and system memory. It is the responsibility of the user to ensure that
+// the appropriate GPU device is selected before each subroutine is called. This
+// is particularly important when using multiple GPU devices on different CPU
+// threads, since active Cuda devices are determined by the cuda runtime on a
+// per-thread basis.
+template <typename T>
+class CudaBuffer {
+ public:
+  CudaBuffer() {}
+  CudaBuffer(const CudaBuffer&) = delete;
+  CudaBuffer& operator=(const CudaBuffer&) = delete;
+
+  ~CudaBuffer() {
+    if (data_ != nullptr) {
+      CHECK_EQ(cudaFree(data_), cudaSuccess);
+    }
+  }
+
+  void Reserve(const size_t size) {
+    if (size > size_) {
+      if (data_ != nullptr) {
+        CHECK_EQ(cudaFree(data_), cudaSuccess);
+      }
+      CHECK_EQ(cudaMalloc(&data_, size * sizeof(T)), cudaSuccess);
+      size_ = size;
+    }
+  }
+
+  void CopyToGpu(const T* data, const size_t size) {
+    Reserve(size);
+    CHECK_EQ(cudaMemcpy(data_, data, size * sizeof(T), cudaMemcpyHostToDevice),
+             cudaSuccess);
+  }
+
+  void CopyToHost(T* data, const size_t size) {
+    CHECK(data_ != nullptr);
+    CHECK_EQ(cudaMemcpy(data, data_, size * sizeof(T), cudaMemcpyDeviceToHost),
+             cudaSuccess);
+  }
+
+  void CopyToGpu(const std::vector<T>& data) {
+    CopyToGpu(data.data(), data.size());
+  }
+
+  T* data() { return data_; }
+  size_t size() const { return size_; }
+
+ private:
+  T* data_ = nullptr;
+  size_t size_ = 0;
+};
+
+#endif  // CERES_NO_CUDA
+
+#endif  // CERES_INTERNAL_CUDA_BUFFER_H_

internal/ceres/cuda_dense_cholesky_test.cc

@@ -0,0 +1,207 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 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: joydeepb@cs.utexas.edu (Joydeep Biswas)
+
+#include <string>
+
+#include "ceres/dense_cholesky.h"
+#include "ceres/internal/eigen.h"
+
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+using std::string;
+
+namespace ceres {
+namespace internal {
+
+#ifndef CERES_NO_CUDA
+
+TEST(DenseCholeskyCudaSolverOld, InvalidOptionOnCreate) {
+  LinearSolver::Options options;
+  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  EXPECT_EQ(dense_cuda_solver, nullptr);
+}
+
+// Tests the CUDA Cholesky solver with a simple 4x4 matrix.
+TEST(DenseCholeskyCudaSolverOld, Cholesky4x4Matrix) {
+  Eigen::Matrix4d A;
+  A <<  4,  12, -16, 0,
+       12,  37, -43, 0,
+      -16, -43,  98, 0,
+        0,   0,   0, 1;
+  const Eigen::Vector4d b = Eigen::Vector4d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
+                                        A.data(),
+                                        &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS);
+  Eigen::Vector4d x = Eigen::Vector4d::Zero();
+  ASSERT_EQ(dense_cuda_solver->Solve(b.data(), x.data(), &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS);
+  EXPECT_NEAR(x(0), 113.75 / 3.0, std::numeric_limits<double>::epsilon() * 10);
+  EXPECT_NEAR(x(1), -31.0 / 3.0, std::numeric_limits<double>::epsilon() * 10);
+  EXPECT_NEAR(x(2), 5.0 / 3.0, std::numeric_limits<double>::epsilon() * 10);
+  EXPECT_NEAR(x(3), 1.0000, std::numeric_limits<double>::epsilon() * 10);
+}
+
+TEST(DenseCholeskyCudaSolverOld, SingularMatrix) {
+  Eigen::Matrix3d A;
+  A <<  1, 0, 0,
+        0, 1, 0,
+        0, 0, 0;
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
+                                        A.data(),
+                                        &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FAILURE);
+}
+
+TEST(DenseCholeskyCudaSolverOld, NegativeMatrix) {
+  Eigen::Matrix3d A;
+  A <<  1, 0, 0,
+        0, 1, 0,
+        0, 0, -1;
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
+                                        A.data(),
+                                        &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FAILURE);
+}
+
+TEST(DenseCholeskyCudaSolverOld, MustFactorizeBeforeSolve) {
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Solve(b.data(), nullptr, &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR);
+}
+
+#ifndef CERES_CUDA_NO_64BIT_SOLVER_API
+// TODO(sameeragarwal): Convert these tests into templated tests so that we use
+// the exact same code for both versions of the test.
+
+TEST(DenseCholeskyCudaSolverNew, InvalidOptionOnCreate) {
+  LinearSolver::Options options;
+  auto dense_cuda_solver = CUDADenseCholesky64Bit::Create(options);
+  EXPECT_EQ(dense_cuda_solver, nullptr);
+}
+
+TEST(DenseCholeskyCudaSolverNew, Cholesky4x4Matrix) {
+  Eigen::Matrix4d A;
+  A <<  4,  12, -16, 0,
+       12,  37, -43, 0,
+      -16, -43,  98, 0,
+        0,   0,   0, 1;
+  const Eigen::Vector4d b = Eigen::Vector4d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky64Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
+                                        A.data(),
+                                        &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS);
+  Eigen::Vector4d x = Eigen::Vector4d::Zero();
+  ASSERT_EQ(dense_cuda_solver->Solve(b.data(), x.data(), &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS);
+  EXPECT_NEAR(x(0), 113.75 / 3.0, std::numeric_limits<double>::epsilon() * 10);
+  EXPECT_NEAR(x(1), -31.0 / 3.0, std::numeric_limits<double>::epsilon() * 10);
+  EXPECT_NEAR(x(2), 5.0 / 3.0, std::numeric_limits<double>::epsilon() * 10);
+  EXPECT_NEAR(x(3), 1.0000, std::numeric_limits<double>::epsilon() * 10);
+}
+
+TEST(DenseCholeskyCudaSolverNew, SingularMatrix) {
+  Eigen::Matrix3d A;
+  A <<  1, 0, 0,
+        0, 1, 0,
+        0, 0, 0;
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky64Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
+                                        A.data(),
+                                        &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FAILURE);
+}
+
+TEST(DenseCholeskyCudaSolverNew, NegativeMatrix) {
+  Eigen::Matrix3d A;
+  A <<  1, 0, 0,
+        0, 1, 0,
+        0, 0, -1;
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky64Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
+                                        A.data(),
+                                        &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FAILURE);
+}
+
+TEST(DenseCholeskyCudaSolverNew, MustFactorizeBeforeSolve) {
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseCholesky64Bit::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  string error_string;
+  ASSERT_EQ(dense_cuda_solver->Solve(b.data(), nullptr, &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR);
+}
+#endif  // CERES_CUDA_NO_64BIT_SOLVER_API
+
+#endif  // CERES_NO_CUDA
+
+}  // namespace internal
+}  // namespace ceres

internal/ceres/dense_cholesky.cc

@@ -33,6 +33,12 @@
 #include <algorithm>
 #include <memory>
 #include <string>
+#include <vector>
+
+#ifndef CERES_NO_CUDA
+#include "cuda_runtime.h"
+#include "cusolverDn.h"
+#endif  // CERES_NO_CUDA
 
 #ifndef CERES_NO_LAPACK
 
@@ -70,6 +76,14 @@ std::unique_ptr<DenseCholesky> DenseCholesky::Create(
       LOG(FATAL) << "Ceres was compiled without support for LAPACK.";
 #endif
 
+    case CUDA:
+#ifndef CERES_NO_CUDA
+      dense_cholesky = CUDADenseCholesky::Create(options);
+      break;
+#else
+      LOG(FATAL) << "Ceres was compiled without support for CUDA.";
+#endif
+
     default:
       LOG(FATAL) << "Unknown dense linear algebra library type : "
                  << DenseLinearAlgebraLibraryTypeToString(
@@ -175,5 +189,362 @@ LinearSolverTerminationType LAPACKDenseCholesky::Solve(const double* rhs,
 
 #endif  // CERES_NO_LAPACK
 
+#ifndef CERES_NO_CUDA
+
+bool CUDADenseCholesky32Bit::Init(std::string* message) {
+  if (cusolverDnCreate(&cusolver_handle_) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnCreate failed.";
+    return false;
+  }
+  if (cudaStreamCreate(&stream_) != cudaSuccess) {
+    *message = "cuSolverDN::cudaStreamCreate failed.";
+    cusolverDnDestroy(cusolver_handle_);
+    return false;
+  }
+  if (cusolverDnSetStream(cusolver_handle_, stream_) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnSetStream failed.";
+    cudaStreamDestroy(stream_);
+    cusolverDnDestroy(cusolver_handle_);
+    return false;
+  }
+  error_.Reserve(1);
+  *message = "CUDADenseCholesky64Bit::Init Success.";
+  return true;
+}
+
+CUDADenseCholesky32Bit::~CUDADenseCholesky32Bit() {
+  if (cusolver_handle_ != nullptr) {
+    CHECK_EQ(cusolverDnDestroy(cusolver_handle_), CUSOLVER_STATUS_SUCCESS);
+    CHECK_EQ(cudaStreamDestroy(stream_), cudaSuccess);
+  }
+}
+
+LinearSolverTerminationType CUDADenseCholesky32Bit::Factorize(
+    int num_cols,
+    double* lhs,
+    std::string* message) {
+  factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  // Allocate GPU memory if necessary.
+  lhs_.Reserve(num_cols * num_cols);
+  num_cols_ = num_cols;
+  // Copy A to GPU.
+  lhs_.CopyToGpu(lhs, num_cols * num_cols);
+  // Allocate scratch space on GPU.
+  int device_workspace_size = 0;
+  if (cusolverDnDpotrf_bufferSize(cusolver_handle_,
+                                      CUBLAS_FILL_MODE_LOWER,
+                                      num_cols,
+                                      lhs_.data(),
+                                      num_cols,
+                                      &device_workspace_size) !=
+      CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDpotrf_bufferSize failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  // Allocate GPU scratch memory.
+  device_workspace_.Reserve(device_workspace_size);
+  // Run the actual factorization (potrf)
+  if (cusolverDnDpotrf(cusolver_handle_,
+                            CUBLAS_FILL_MODE_LOWER,
+                            num_cols,
+                            lhs_.data(),
+                            num_cols,
+                            reinterpret_cast<double*>(device_workspace_.data()),
+                            device_workspace_.size(),
+                            error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDpotrf failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_), cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  int error = 0;
+  // Check for errors.
+  error_.CopyToHost(&error, 1);
+  if (error < 0) {
+    LOG(FATAL) << "Congratulations, you found a bug in Ceres - "
+               << "please report it. "
+               << "cuSolverDN::cusolverDnXpotrf fatal error."
+               << "Argument: " << -error << " is invalid.";
+    // The following line is unreachable, but return failure just to be
+    // pedantic, since the compiler does not know that.
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  } else if (error > 0) {
+    *message = StringPrintf(
+        "cuSolverDN::cusolverDnDpotrf numerical failure. "
+        "The leading minor of order %d is not positive definite.",
+        error);
+    factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_FAILURE;
+    return LinearSolverTerminationType::LINEAR_SOLVER_FAILURE;
+  }
+  *message = "Success";
+  factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+  return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+}
+
+LinearSolverTerminationType CUDADenseCholesky32Bit::Solve(
+    const double* rhs,
+    double* solution,
+    std::string* message) {
+  if (factorize_result_ != LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS) {
+    *message = "Factorize did not complete succesfully previously.";
+    return factorize_result_;
+  }
+  // Copy RHS to GPU.
+  rhs_.CopyToGpu(rhs, num_cols_);
+  // Solve the system.
+  if (cusolverDnDpotrs(cusolver_handle_,
+                            CUBLAS_FILL_MODE_LOWER,
+                            num_cols_,
+                            1,
+                            lhs_.data(),
+                            num_cols_,
+                            rhs_.data(),
+                            num_cols_,
+                            error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDpotrs failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_), cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  // Check for errors.
+  int error = 0;
+  // Copy X from GPU to host.
+  rhs_.CopyToHost(solution, num_cols_);
+  if (error != 0) {
+    LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+               << "Please report it."
+               << "cuSolverDN::cusolverDnDpotrs fatal error."
+               << "Argument: " << -error << " is invalid.";
+  }
+  // Copy error variable from GPU to host.
+  error_.CopyToHost(&error, 1);
+  *message = "Success";
+  return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+}
+
+std::unique_ptr<CUDADenseCholesky32Bit> CUDADenseCholesky32Bit::Create(
+  const LinearSolver::Options& options) {
+  if (options.dense_linear_algebra_library_type != CUDA) {
+    // The user called the wrong factory method.
+    return nullptr;
+  }
+  CUDADenseCholesky32Bit* cuda_dense_cholesky = new CUDADenseCholesky32Bit();
+  std::string cuda_error;
+  if (cuda_dense_cholesky->Init(&cuda_error)) {
+    return std::unique_ptr<CUDADenseCholesky32Bit>(cuda_dense_cholesky);
+  }
+  // Initialization failed, destroy the object and return a nullptr.
+  delete cuda_dense_cholesky;
+  LOG(ERROR) << "CUDADenseCholesky32Bit::Init failed: " << cuda_error;
+  return std::unique_ptr<CUDADenseCholesky32Bit>(nullptr);
+}
+
+std::unique_ptr<CUDADenseCholesky64Bit> CUDADenseCholesky64Bit::Create(
+  const LinearSolver::Options& options) {
+  if (options.dense_linear_algebra_library_type != CUDA) {
+    // The user called the wrong factory method.
+    return nullptr;
+  }
+  CUDADenseCholesky64Bit* cuda_dense_cholesky = new CUDADenseCholesky64Bit();
+  std::string cuda_error;
+  if (cuda_dense_cholesky->Init(&cuda_error)) {
+    return std::unique_ptr<CUDADenseCholesky64Bit>(cuda_dense_cholesky);
+  }
+  // Initialization failed, destroy the object and return a nullptr.
+  delete cuda_dense_cholesky;
+  LOG(ERROR) << "CUDADenseCholesky64Bit::Init failed: " << cuda_error;
+  return std::unique_ptr<CUDADenseCholesky64Bit>(nullptr);
+}
+
+#ifdef CERES_CUDA_NO_64BIT_SOLVER_API
+// CUDA < 11.1 did not have the 64-bit APIs, so the implementation of the new
+// interface will just generate a fatal error.
+
+CUDADenseCholesky64Bit::~CUDADenseCholesky64Bit() {}
+
+bool CUDADenseCholesky64Bit::Init(std::string* message) {
+  *message = "Cannot use CUDADenseCholesky64Bit with CUDA < 11.1.";
+  return false;
+}
+
+LinearSolverTerminationType CUDADenseCholesky64Bit::Factorize(
+    int,
+    double*,
+    std::string*) {
+  return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+}
+
+LinearSolverTerminationType CUDADenseCholesky64Bit::Solve(
+    const double*,
+    double*,
+    std::string*) {
+  return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+}
+
+#else  // CERES_CUDA_NO_64BIT_SOLVER_API
+
+bool CUDADenseCholesky64Bit::Init(std::string* message) {
+  if (cusolverDnCreate(&cusolver_handle_) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnCreate failed.";
+    return false;
+  }
+  if (cudaStreamCreate(&stream_) != cudaSuccess) {
+    *message = "cuSolverDN::cudaStreamCreate failed.";
+    cusolverDnDestroy(cusolver_handle_);
+    return false;
+  }
+  if (cusolverDnSetStream(cusolver_handle_, stream_) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnSetStream failed.";
+    cudaStreamDestroy(stream_);
+    cusolverDnDestroy(cusolver_handle_);
+    return false;
+  }
+  error_.Reserve(1);
+  *message = "CUDADenseCholesky64Bit::Init Success.";
+  return true;
+}
+
+CUDADenseCholesky64Bit::~CUDADenseCholesky64Bit() {
+  if (cusolver_handle_ != nullptr) {
+    CHECK_EQ(cusolverDnDestroy(cusolver_handle_), CUSOLVER_STATUS_SUCCESS);
+    CHECK_EQ(cudaStreamDestroy(stream_), cudaSuccess);
+  }
+}
+
+LinearSolverTerminationType CUDADenseCholesky64Bit::Factorize(
+    int num_cols,
+    double* lhs,
+    std::string* message) {
+  factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  // Allocate GPU memory if necessary.
+  lhs_.Reserve(num_cols * num_cols);
+  num_cols_ = num_cols;
+  // Copy A to GPU.
+  lhs_.CopyToGpu(lhs, num_cols * num_cols);
+  // Allocate scratch space on GPU.
+  size_t host_workspace_size = 0;
+  size_t device_workspace_size = 0;
+  if (cusolverDnXpotrf_bufferSize(cusolver_handle_,
+                                      nullptr,
+                                      CUBLAS_FILL_MODE_LOWER,
+                                      num_cols,
+                                      CUDA_R_64F,
+                                      lhs_.data(),
+                                      num_cols,
+                                      CUDA_R_64F,
+                                      &device_workspace_size,
+                                      &host_workspace_size) !=
+      CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnXpotrf_bufferSize failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  // Allocate host scratch memory.
+  host_workspace_.resize(host_workspace_size);
+  // Allocate GPU scratch memory.
+  device_workspace_.Reserve(device_workspace_size);
+  // Run the actual factorization (potrf)
+  if (cusolverDnXpotrf(cusolver_handle_,
+                            nullptr,
+                            CUBLAS_FILL_MODE_LOWER,
+                            num_cols,
+                            CUDA_R_64F,
+                            lhs_.data(),
+                            num_cols,
+                            CUDA_R_64F,
+                            device_workspace_.data(),
+                            device_workspace_.size(),
+                            host_workspace_.data(),
+                            host_workspace_.size(),
+                            error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnXpotrf failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_) != cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  int error = 0;
+  // Check for errors.
+  error_.CopyToHost(&error, 1);
+  if (error < 0) {
+    LOG(FATAL) << "Congratulations, you found a bug in Ceres - "
+               << "please report it. "
+               << "cuSolverDN::cusolverDnXpotrf fatal error."
+               << "Argument: " << -error << " is invalid.";
+    // The following line is unreachable, but return failure just to be
+    // pedantic, since the compiler does not know that.
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  } else if (error > 0) {
+    *message = StringPrintf(
+        "cuSolverDN::cusolverDnXpotrf numerical failure. "
+        "The leading minor of order %d is not positive definite.",
+        error);
+    factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_FAILURE;
+    return LinearSolverTerminationType::LINEAR_SOLVER_FAILURE;
+  }
+
+  *message = "Success";
+  factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+  return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+}
+
+LinearSolverTerminationType CUDADenseCholesky64Bit::Solve(
+    const double* rhs,
+    double* solution,
+    std::string* message) {
+  if (factorize_result_ != LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS) {
+    *message = "Factorize did not complete succesfully previously.";
+    return factorize_result_;
+  }
+  // Copy RHS to GPU.
+  rhs_.CopyToGpu(rhs, num_cols_);
+  // Solve the system.
+  if (cusolverDnXpotrs(cusolver_handle_,
+                            nullptr,
+                            CUBLAS_FILL_MODE_LOWER,
+                            num_cols_,
+                            1,
+                            CUDA_R_64F,
+                            lhs_.data(),
+                            num_cols_,
+                            CUDA_R_64F,
+                            rhs_.data(),
+                            num_cols_,
+                            error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnXpotrs failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_) != cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  // Check for errors.
+  int error = 0;
+  // Copy error variable from GPU to host.
+  error_.CopyToHost(&error, 1);
+  if (error != 0) {
+    LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+               << "Please report it."
+               << "cuSolverDN::cusolverDnXpotrs fatal error."
+               << "Argument: " << -error << " is invalid.";
+  }
+  // Copy X from GPU to host.
+  rhs_.CopyToHost(solution, num_cols_);
+  *message = "Success";
+  return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+}
+
+#endif // CERES_CUDA_NO_64BIT_SOLVER_API
+
+#endif  // CERES_NO_CUDA
+
 }  // namespace internal
 }  // namespace ceres

internal/ceres/dense_cholesky.h

@@ -37,10 +37,16 @@
 // clang-format on
 
 #include <memory>
+#include <vector>
 
 #include "Eigen/Dense"
+#include "ceres/cuda_buffer.h"
 #include "ceres/linear_solver.h"
 #include "glog/logging.h"
+#ifndef CERES_NO_CUDA
+#include "cuda_runtime.h"
+#include "cusolverDn.h"
+#endif  // CERES_NO_CUDA
 
 namespace ceres {
 namespace internal {
@@ -129,6 +135,107 @@ class CERES_EXPORT_INTERNAL LAPACKDenseCholesky : public DenseCholesky {
 };
 #endif  // CERES_NO_LAPACK
 
+#ifndef CERES_NO_CUDA
+// Implementation of DenseCholesky using the cuSolver library v.11.0 or older,
+// using the legacy cuSolverDn interface.
+class CERES_EXPORT_INTERNAL CUDADenseCholesky32Bit : public DenseCholesky {
+ public:
+  static std::unique_ptr<CUDADenseCholesky32Bit> Create(
+      const LinearSolver::Options& options);
+  ~CUDADenseCholesky32Bit() override;
+  CUDADenseCholesky32Bit(const CUDADenseCholesky32Bit&) = delete;
+  CUDADenseCholesky32Bit& operator=(const CUDADenseCholesky32Bit&) = delete;
+  LinearSolverTerminationType Factorize(int num_cols,
+                                        double* lhs,
+                                        std::string* message) override;
+  LinearSolverTerminationType Solve(const double* rhs,
+                                    double* solution,
+                                    std::string* message) override;
+
+ private:
+  CUDADenseCholesky32Bit() {}
+  // Initializes the cuSolverDN context, creates an asynchronous stream, and
+  // associates the stream with cuSolverDN. Returns true iff initialization was
+  // successful, else it returns false and a human-readable error message is
+  // returned.
+  bool Init(std::string* message);
+
+  // Handle to the cuSOLVER context.
+  cusolverDnHandle_t cusolver_handle_ = nullptr;
+  // CUDA device stream.
+  cudaStream_t stream_ = nullptr;
+  // Number of columns in the A matrix, to be cached between calls to *Factorize
+  // and *Solve.
+  size_t num_cols_ = 0;
+  // GPU memory allocated for the A matrix (lhs matrix).
+  CudaBuffer<double> lhs_;
+  // GPU memory allocated for the B matrix (rhs vector).
+  CudaBuffer<double> rhs_;
+  // Scratch space for cuSOLVER on the GPU.
+  CudaBuffer<uint8_t> device_workspace_;
+  // Required for error handling with cuSOLVER.
+  CudaBuffer<int> error_;
+  // Cache the result of Factorize to ensure that when Solve is called, the
+  // factiorization of lhs is valid.
+  LinearSolverTerminationType factorize_result_ =
+      LINEAR_SOLVER_FATAL_ERROR;
+};
+
+// Implementation of DenseCholesky using the cuSolver library v.11.1 or newer,
+// using the 64-bit cuSolverDn interface.
+class CERES_EXPORT_INTERNAL CUDADenseCholesky64Bit : public DenseCholesky {
+ public:
+  static std::unique_ptr<CUDADenseCholesky64Bit> Create(
+      const LinearSolver::Options& options);
+  ~CUDADenseCholesky64Bit() override;
+  CUDADenseCholesky64Bit(const CUDADenseCholesky64Bit&) = delete;
+  CUDADenseCholesky64Bit& operator=(const CUDADenseCholesky64Bit&) = delete;
+  LinearSolverTerminationType Factorize(int num_cols,
+                                        double* lhs,
+                                        std::string* message) override;
+  LinearSolverTerminationType Solve(const double* rhs,
+                                    double* solution,
+                                    std::string* message) override;
+
+ private:
+  CUDADenseCholesky64Bit() {}
+  // Initializes the cuSolverDN context, creates an asynchronous stream, and
+  // associates the stream with cuSolverDN. Returns true iff initialization was
+  // successful, else it returns false and a human-readable error message is
+  // returned.
+  bool Init(std::string* message);
+
+  // Handle to the cuSOLVER context.
+  cusolverDnHandle_t cusolver_handle_ = nullptr;
+  // CUDA device stream.
+  cudaStream_t stream_ = nullptr;
+  // Number of columns in the A matrix, to be cached between calls to Factorize
+  // and Solve.
+  size_t num_cols_ = 0;
+  // GPU memory allocated for the A matrix (lhs matrix).
+  CudaBuffer<double> lhs_;
+  // GPU memory allocated for the B matrix (rhs vector).
+  CudaBuffer<double> rhs_;
+  // Workspace for cuSOLVER on the GPU.
+  CudaBuffer<uint8_t> device_workspace_;
+  // Workspace for cuSOLVER on the host.
+  std::vector<double> host_workspace_;
+  // Required for error handling with cuSOLVER.
+  CudaBuffer<int> error_;
+  // Cache the result of Factorize to ensure that when Solve is called, the
+  // factiorization of lhs is valid.
+  LinearSolverTerminationType factorize_result_ =
+      LINEAR_SOLVER_FATAL_ERROR;
+};
+
+#ifdef CERES_CUDA_NO_64BIT_SOLVER_API
+using CUDADenseCholesky = CUDADenseCholesky32Bit;
+#else
+using CUDADenseCholesky = CUDADenseCholesky64Bit;
+#endif
+
+#endif  // CERES_NO_CUDA
+
 }  // namespace internal
 }  // namespace ceres
 

internal/ceres/dense_cholesky_test.cc

@@ -92,17 +92,20 @@ TEST_P(DenseCholeskyTest, FactorAndSolve) {
   }
 }
 
+INSTANTIATE_TEST_SUITE_P(
+    _,
+    DenseCholeskyTest,
+    ::testing::Values(
+        EIGEN,
 #ifndef CERES_NO_LAPACK
-INSTANTIATE_TEST_SUITE_P(_,
-                         DenseCholeskyTest,
-                         ::testing::Values(EIGEN, LAPACK),
-                         ParamInfoToString);
-#else
-INSTANTIATE_TEST_SUITE_P(_,
-                         DenseCholeskyTest,
-                         ::testing::Values(EIGEN),
-                         ParamInfoToString);
+        LAPACK,
 #endif
+#ifndef CERES_NO_CUDA
+        CUDA
+#endif
+    ),
+    ParamInfoToString);
+
 
 }  // namespace internal
 }  // namespace ceres

internal/ceres/dense_linear_solver_benchmark.cc

@@ -95,7 +95,12 @@ BENCHMARK_TEMPLATE2(BM_DenseSolver, ceres::LAPACK, ceres::DENSE_QR)
     ->Apply(MatrixSizes);
 BENCHMARK_TEMPLATE2(BM_DenseSolver, ceres::LAPACK, ceres::DENSE_NORMAL_CHOLESKY)
     ->Apply(MatrixSizes);
-#endif
+#endif  // CERES_NO_LAPACK
+
+#ifndef CERES_NO_CUDA
+BENCHMARK(BM_DenseSolver<ceres::CUDA, ceres::DENSE_NORMAL_CHOLESKY>)
+    ->Apply(MatrixSizes);
+#endif  // CERES_NO_CUDA
 
 }  // namespace internal
 }  // namespace ceres

internal/ceres/solver_utils.cc

@@ -36,6 +36,9 @@
 #include "ceres/internal/config.h"
 #include "ceres/internal/port.h"
 #include "ceres/version.h"
+#ifndef CERES_NO_CUDA
+#include "cuda_runtime.h"
+#endif  // CERES_NO_CUDA
 
 namespace ceres {
 namespace internal {
@@ -87,6 +90,10 @@ std::string VersionString() {
   value += "-no_custom_blas";
 #endif
 
+#ifndef CERES_NO_CUDA
+  value += "-cuda-(" + std::to_string(CUDART_VERSION) + ")";
+#endif
+
   return value;
 }
 

internal/ceres/types.cc

@@ -128,6 +128,7 @@ const char* DenseLinearAlgebraLibraryTypeToString(
   switch (type) {
     CASESTR(EIGEN);
     CASESTR(LAPACK);
+    CASESTR(CUDA);
     default:
       return "UNKNOWN";
   }
@@ -138,6 +139,7 @@ bool StringToDenseLinearAlgebraLibraryType(
   UpperCase(&value);
   STRENUM(EIGEN);
   STRENUM(LAPACK);
+  STRENUM(CUDA);
   return false;
 }