diff --git a/CMakeLists.txt b/CMakeLists.txt
index 48ba201..69432f7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -270,14 +270,6 @@
         "${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)
diff --git a/examples/nist.cc b/examples/nist.cc
index a3430ca..cb0d1aa0 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -101,7 +101,7 @@
               "and cgnr");
 DEFINE_string(dense_linear_algebra_library,
               "eigen",
-              "Options are: eigen and lapack.");
+              "Options are: eigen, lapack, and cuda.");
 DEFINE_string(preconditioner, "jacobi", "Options are: identity, jacobi");
 DEFINE_string(line_search,
               "wolfe",
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 166eb19..6c6c3c1 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -143,6 +143,7 @@
   list(APPEND
        CERES_LIBRARY_PRIVATE_DEPENDENCIES
        ${CUDA_LIBRARIES}
+       ${CUDA_cublas_LIBRARY}
        ${CUDA_cusolver_LIBRARY}
        ${CUDA_cusparse_LIBRARY})
 endif (CUDA_FOUND)
@@ -479,6 +480,7 @@
   ceres_test(covariance)
   ceres_test(cubic_interpolation)
   ceres_test(cuda_dense_cholesky)
+  ceres_test(cuda_dense_qr)
   ceres_test(dense_linear_solver)
   ceres_test(dense_cholesky)
   ceres_test(dense_qr)
diff --git a/internal/ceres/cuda_buffer.h b/internal/ceres/cuda_buffer.h
index 44e7de3..61f9752 100644
--- a/internal/ceres/cuda_buffer.h
+++ b/internal/ceres/cuda_buffer.h
@@ -45,7 +45,8 @@
 // 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.
+// per-thread basis. Note that unless otherwise specified, all methods use the
+// default stream, and are synchronous.
 template <typename T>
 class CudaBuffer {
  public:
@@ -59,6 +60,8 @@
     }
   }
 
+  // Grow the GPU memory buffer if needed to accommodate data of the specified
+  // size
   void Reserve(const size_t size) {
     if (size > size_) {
       if (data_ != nullptr) {
@@ -69,12 +72,21 @@
     }
   }
 
-  void CopyToGpu(const T* data, const size_t size) {
+  // Perform an asynchronous copy from CPU memory to GPU memory using the stream
+  // provided.
+  void CopyToGpuAsync(const T* data, const size_t size, cudaStream_t stream) {
     Reserve(size);
-    CHECK_EQ(cudaMemcpy(data_, data, size * sizeof(T), cudaMemcpyHostToDevice),
+    CHECK_EQ(cudaMemcpyAsync(data_,
+                             data,
+                             size * sizeof(T),
+                             cudaMemcpyHostToDevice,
+                             stream),
              cudaSuccess);
   }
 
+  // Copy data from the GPU to CPU memory. This is necessarily synchronous since
+  // any potential GPU kernels that may be writing to the buffer must finish
+  // before the transfer happens.
   void CopyToHost(T* data, const size_t size) {
     CHECK(data_ != nullptr);
     CHECK_EQ(cudaMemcpy(data, data_, size * sizeof(T), cudaMemcpyDeviceToHost),
diff --git a/internal/ceres/cuda_dense_cholesky_test.cc b/internal/ceres/cuda_dense_cholesky_test.cc
index 7d39b0e..6c4dcc3 100644
--- a/internal/ceres/cuda_dense_cholesky_test.cc
+++ b/internal/ceres/cuda_dense_cholesky_test.cc
@@ -36,21 +36,19 @@
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
-using std::string;
-
 namespace ceres {
 namespace internal {
 
 #ifndef CERES_NO_CUDA
 
-TEST(DenseCholeskyCudaSolverOld, InvalidOptionOnCreate) {
+TEST(CUDADenseCholesky, InvalidOptionOnCreate) {
   LinearSolver::Options options;
-  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  auto dense_cuda_solver = CUDADenseCholesky::Create(options);
   EXPECT_EQ(dense_cuda_solver, nullptr);
 }
 
 // Tests the CUDA Cholesky solver with a simple 4x4 matrix.
-TEST(DenseCholeskyCudaSolverOld, Cholesky4x4Matrix) {
+TEST(CUDADenseCholesky, Cholesky4x4Matrix) {
   Eigen::Matrix4d A;
   A <<  4,  12, -16, 0,
        12,  37, -43, 0,
@@ -59,9 +57,9 @@
   const Eigen::Vector4d b = Eigen::Vector4d::Ones();
   LinearSolver::Options options;
   options.dense_linear_algebra_library_type = CUDA;
-  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  auto dense_cuda_solver = CUDADenseCholesky::Create(options);
   ASSERT_NE(dense_cuda_solver, nullptr);
-  string error_string;
+  std::string error_string;
   ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
                                         A.data(),
                                         &error_string),
@@ -75,7 +73,7 @@
   EXPECT_NEAR(x(3), 1.0000, std::numeric_limits<double>::epsilon() * 10);
 }
 
-TEST(DenseCholeskyCudaSolverOld, SingularMatrix) {
+TEST(CUDADenseCholesky, SingularMatrix) {
   Eigen::Matrix3d A;
   A <<  1, 0, 0,
         0, 1, 0,
@@ -83,16 +81,16 @@
   const Eigen::Vector3d b = Eigen::Vector3d::Ones();
   LinearSolver::Options options;
   options.dense_linear_algebra_library_type = CUDA;
-  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  auto dense_cuda_solver = CUDADenseCholesky::Create(options);
   ASSERT_NE(dense_cuda_solver, nullptr);
-  string error_string;
+  std::string error_string;
   ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
                                         A.data(),
                                         &error_string),
             LinearSolverTerminationType::LINEAR_SOLVER_FAILURE);
 }
 
-TEST(DenseCholeskyCudaSolverOld, NegativeMatrix) {
+TEST(CUDADenseCholesky, NegativeMatrix) {
   Eigen::Matrix3d A;
   A <<  1, 0, 0,
         0, 1, 0,
@@ -100,107 +98,26 @@
   const Eigen::Vector3d b = Eigen::Vector3d::Ones();
   LinearSolver::Options options;
   options.dense_linear_algebra_library_type = CUDA;
-  auto dense_cuda_solver = CUDADenseCholesky32Bit::Create(options);
+  auto dense_cuda_solver = CUDADenseCholesky::Create(options);
   ASSERT_NE(dense_cuda_solver, nullptr);
-  string error_string;
+  std::string error_string;
   ASSERT_EQ(dense_cuda_solver->Factorize(A.cols(),
                                         A.data(),
                                         &error_string),
             LinearSolverTerminationType::LINEAR_SOLVER_FAILURE);
 }
 
-TEST(DenseCholeskyCudaSolverOld, MustFactorizeBeforeSolve) {
+TEST(CUDADenseCholesky, 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);
+  auto dense_cuda_solver = CUDADenseCholesky::Create(options);
   ASSERT_NE(dense_cuda_solver, nullptr);
-  string error_string;
+  std::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
diff --git a/internal/ceres/cuda_dense_qr_test.cc b/internal/ceres/cuda_dense_qr_test.cc
new file mode 100644
index 0000000..15ba00a
--- /dev/null
+++ b/internal/ceres/cuda_dense_qr_test.cc
@@ -0,0 +1,164 @@
+// 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_qr.h"
+#include "ceres/internal/eigen.h"
+
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+#ifndef CERES_NO_CUDA
+
+TEST(CUDADenseQR, InvalidOptionOnCreate) {
+  LinearSolver::Options options;
+  auto dense_cuda_solver = CUDADenseQR::Create(options);
+  EXPECT_EQ(dense_cuda_solver, nullptr);
+}
+
+// Tests the CUDA QR solver with a simple 4x4 matrix.
+TEST(CUDADenseQR, QR4x4Matrix) {
+  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 = CUDADenseQR::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  std::string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.rows(),
+                                         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);
+  // Empirically observed accuracy of cuSolverDN's QR solver.
+  const double kEpsilon = 1e-11;
+  EXPECT_NEAR(x(0), 113.75 / 3.0, kEpsilon);
+  EXPECT_NEAR(x(1), -31.0 / 3.0, kEpsilon);
+  EXPECT_NEAR(x(2), 5.0 / 3.0, kEpsilon);
+  EXPECT_NEAR(x(3), 1.0000, kEpsilon);
+}
+
+// Tests the CUDA QR solver with a simple 4x4 matrix.
+TEST(CUDADenseQR, QR4x2Matrix) {
+  Eigen::Matrix<double, 4, 2> A;
+  A <<  4,  12,
+       12,  37,
+      -16, -43,
+        0,   0;
+  const std::vector<double> b(4, 1.0);
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseQR::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  std::string error_string;
+  ASSERT_EQ(dense_cuda_solver->Factorize(A.rows(),
+                                         A.cols(),
+                                         A.data(),
+                                         &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS);
+  std::vector<double> x(2, 0);
+  ASSERT_EQ(dense_cuda_solver->Solve(b.data(), x.data(), &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS);
+  // Empirically observed accuracy of cuSolverDN's QR solver.
+  const double kEpsilon = 1e-11;
+  // Solution values computed with Octave.
+  EXPECT_NEAR(x[0], -1.143410852713177, kEpsilon);
+  EXPECT_NEAR(x[1], 0.4031007751937981, kEpsilon);
+}
+
+TEST(CUDADenseQR, MustFactorizeBeforeSolve) {
+  const Eigen::Vector3d b = Eigen::Vector3d::Ones();
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = CUDA;
+  auto dense_cuda_solver = CUDADenseQR::Create(options);
+  ASSERT_NE(dense_cuda_solver, nullptr);
+  std::string error_string;
+  ASSERT_EQ(dense_cuda_solver->Solve(b.data(), nullptr, &error_string),
+            LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR);
+}
+
+TEST(CUDADenseQR, Randomized1600x100Tests) {
+  const int kNumRows = 1600;
+  const int kNumCols = 100;
+  using LhsType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
+  using RhsType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+  using SolutionType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+
+  LinearSolver::Options options;
+  options.dense_linear_algebra_library_type = ceres::CUDA;
+  std::unique_ptr<DenseQR> dense_qr = CUDADenseQR::Create(options);
+
+  const int kNumTrials = 100;
+  const int kMinNumCols = 1;
+  const int kMaxNumCols = 10;
+  const int kMinRowsFactor = 1;
+  const int kMaxRowsFactor = 3;
+  for (int i = 0; i < kNumTrials; ++i) {
+    LhsType lhs = LhsType::Random(kNumRows, kNumCols);
+    SolutionType x_expected = SolutionType::Random(kNumCols);
+    RhsType rhs = lhs * x_expected;
+    SolutionType x_computed = SolutionType::Zero(kNumCols);
+    // Sanity check the random matrix sizes.
+    EXPECT_EQ(lhs.rows(), kNumRows);
+    EXPECT_EQ(lhs.cols(), kNumCols);
+    EXPECT_EQ(rhs.rows(), kNumRows);
+    EXPECT_EQ(rhs.cols(), 1);
+    EXPECT_EQ(x_expected.rows(), kNumCols);
+    EXPECT_EQ(x_expected.cols(), 1);
+    EXPECT_EQ(x_computed.rows(), kNumCols);
+    EXPECT_EQ(x_computed.cols(), 1);
+    LinearSolver::Summary summary;
+    summary.termination_type = dense_qr->FactorAndSolve(kNumRows,
+                                                        kNumCols,
+                                                        lhs.data(),
+                                                        rhs.data(),
+                                                        x_computed.data(),
+                                                        &summary.message);
+    ASSERT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
+    ASSERT_NEAR((x_computed - x_expected).norm() / x_expected.norm(),
+                0.0,
+                std::numeric_limits<double>::epsilon() * 400);
+  }
+}
+#endif  // CERES_NO_CUDA
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/dense_cholesky.cc b/internal/ceres/dense_cholesky.cc
index af555e7..2d8c2da 100644
--- a/internal/ceres/dense_cholesky.cc
+++ b/internal/ceres/dense_cholesky.cc
@@ -193,13 +193,14 @@
 
 #ifndef CERES_NO_CUDA
 
-bool CUDADenseCholesky32Bit::Init(std::string* message) {
+bool CUDADenseCholesky::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.";
+  if (cudaStreamCreateWithFlags(&stream_, cudaStreamNonBlocking) !=
+      cudaSuccess) {
+    *message = "CUDA::cudaStreamCreateWithFlags failed.";
     cusolverDnDestroy(cusolver_handle_);
     return false;
   }
@@ -211,26 +212,23 @@
     return false;
   }
   error_.Reserve(1);
-  *message = "CUDADenseCholesky64Bit::Init Success.";
+  *message = "CUDADenseCholesky::Init Success.";
   return true;
 }
 
-CUDADenseCholesky32Bit::~CUDADenseCholesky32Bit() {
+CUDADenseCholesky::~CUDADenseCholesky() {
   if (cusolver_handle_ != nullptr) {
     CHECK_EQ(cusolverDnDestroy(cusolver_handle_), CUSOLVER_STATUS_SUCCESS);
     CHECK_EQ(cudaStreamDestroy(stream_), cudaSuccess);
   }
 }
 
-LinearSolverTerminationType CUDADenseCholesky32Bit::Factorize(
+LinearSolverTerminationType CUDADenseCholesky::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.
+  lhs_.CopyToGpuAsync(lhs, num_cols * num_cols, stream_);
   int device_workspace_size = 0;
   if (cusolverDnDpotrf_bufferSize(cusolver_handle_,
                                   CUBLAS_FILL_MODE_LOWER,
@@ -242,9 +240,7 @@
     *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,
@@ -256,13 +252,12 @@
     *message = "cuSolverDN::cusolverDnDpotrf failed.";
     return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
   }
-  if (cudaDeviceSynchronize() != cudaSuccess || cudaStreamSynchronize(stream_),
-      cudaSuccess) {
+  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 - "
@@ -285,15 +280,13 @@
   return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
 }
 
-LinearSolverTerminationType CUDADenseCholesky32Bit::Solve(
+LinearSolverTerminationType CUDADenseCholesky::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.
+  rhs_.CopyToGpuAsync(rhs, num_cols_, stream_);
   if (cusolverDnDpotrs(cusolver_handle_,
                        CUBLAS_FILL_MODE_LOWER,
                        num_cols_,
@@ -306,237 +299,41 @@
     *message = "cuSolverDN::cusolverDnDpotrs failed.";
     return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
   }
-  if (cudaDeviceSynchronize() != cudaSuccess || cudaStreamSynchronize(stream_),
-      cudaSuccess) {
+  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_);
+  error_.CopyToHost(&error, 1);
   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
+std::unique_ptr<CUDADenseCholesky> CUDADenseCholesky::Create(
+    const LinearSolver::Options& options) {
+  if (options.dense_linear_algebra_library_type != CUDA) {
+    // The user called the wrong factory method.
+    return nullptr;
+  }
+  auto cuda_dense_cholesky =
+      std::unique_ptr<CUDADenseCholesky>(new CUDADenseCholesky());
+  std::string cuda_error;
+  if (cuda_dense_cholesky->Init(&cuda_error)) {
+    return cuda_dense_cholesky;
+  }
+  // Initialization failed, destroy the object (done automatically) and return a
+  // nullptr.
+  LOG(ERROR) << "CUDADenseCholesky::Init failed: " << cuda_error;
+  return nullptr;
+}
 
 #endif  // CERES_NO_CUDA
 
diff --git a/internal/ceres/dense_cholesky.h b/internal/ceres/dense_cholesky.h
index f7e48d4..49d780c 100644
--- a/internal/ceres/dense_cholesky.h
+++ b/internal/ceres/dense_cholesky.h
@@ -134,15 +134,15 @@
 #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_NO_EXPORT CUDADenseCholesky32Bit : public DenseCholesky {
+// CUDA implementation of DenseCholesky using the cuSolverDN library using the
+// 32-bit legacy interface for maximum compatibility.
+class CERES_NO_EXPORT CUDADenseCholesky : public DenseCholesky {
  public:
-  static std::unique_ptr<CUDADenseCholesky32Bit> Create(
+  static std::unique_ptr<CUDADenseCholesky> Create(
       const LinearSolver::Options& options);
-  ~CUDADenseCholesky32Bit() override;
-  CUDADenseCholesky32Bit(const CUDADenseCholesky32Bit&) = delete;
-  CUDADenseCholesky32Bit& operator=(const CUDADenseCholesky32Bit&) = delete;
+  ~CUDADenseCholesky() override;
+  CUDADenseCholesky(const CUDADenseCholesky&) = delete;
+  CUDADenseCholesky& operator=(const CUDADenseCholesky&) = delete;
   LinearSolverTerminationType Factorize(int num_cols,
                                         double* lhs,
                                         std::string* message) override;
@@ -151,7 +151,7 @@
                                     std::string* message) override;
 
  private:
-  CUDADenseCholesky32Bit() = default;
+  CUDADenseCholesky() = default;
   // 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
@@ -179,59 +179,6 @@
       LINEAR_SOLVER_FATAL_ERROR;
 };
 
-// Implementation of DenseCholesky using the cuSolver library v.11.1 or newer,
-// using the 64-bit cuSolverDn interface.
-class CERES_NO_EXPORT 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() = default;
-  // 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
diff --git a/internal/ceres/dense_linear_solver_benchmark.cc b/internal/ceres/dense_linear_solver_benchmark.cc
index 2af9fa8..ba171bf 100644
--- a/internal/ceres/dense_linear_solver_benchmark.cc
+++ b/internal/ceres/dense_linear_solver_benchmark.cc
@@ -98,7 +98,9 @@
 #endif  // CERES_NO_LAPACK
 
 #ifndef CERES_NO_CUDA
-BENCHMARK(BM_DenseSolver<ceres::CUDA, ceres::DENSE_NORMAL_CHOLESKY>)
+BENCHMARK_TEMPLATE2(BM_DenseSolver, ceres::CUDA, ceres::DENSE_NORMAL_CHOLESKY)
+    ->Apply(MatrixSizes);
+BENCHMARK_TEMPLATE2(BM_DenseSolver, ceres::CUDA, ceres::DENSE_QR)
     ->Apply(MatrixSizes);
 #endif  // CERES_NO_CUDA
 
diff --git a/internal/ceres/dense_qr.cc b/internal/ceres/dense_qr.cc
index 33a7cbd..77f04ad 100644
--- a/internal/ceres/dense_qr.cc
+++ b/internal/ceres/dense_qr.cc
@@ -33,6 +33,10 @@
 #include <algorithm>
 #include <memory>
 #include <string>
+#ifndef CERES_NO_CUDA
+#include "cusolverDn.h"
+#include "cublas_v2.h"
+#endif  // CERES_NO_CUDA
 
 #ifndef CERES_NO_LAPACK
 
@@ -124,6 +128,14 @@
       LOG(FATAL) << "Ceres was compiled without support for LAPACK.";
 #endif
 
+    case CUDA:
+#ifndef CERES_NO_CUDA
+      dense_qr = CUDADenseQR::Create(options);
+      break;
+#else
+      LOG(FATAL) << "Ceres was compiled without support for CUDA.";
+#endif
+
     default:
       LOG(FATAL) << "Unknown dense linear algebra library type : "
                  << DenseLinearAlgebraLibraryTypeToString(
@@ -296,5 +308,207 @@
 
 #endif  // CERES_NO_LAPACK
 
+#ifndef CERES_NO_CUDA
+
+bool CUDADenseQR::Init(std::string* message) {
+  if (cublasCreate(&cublas_handle_) != CUBLAS_STATUS_SUCCESS) {
+    *message = "cuBLAS::cublasCreate failed.";
+    return false;
+  }
+  if (cusolverDnCreate(&cusolver_handle_) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnCreate failed.";
+    return false;
+  }
+  if (cudaStreamCreateWithFlags(&stream_, cudaStreamNonBlocking) !=
+      cudaSuccess) {
+    *message = "CUDA::cudaStreamCreateWithFlags failed.";
+    cusolverDnDestroy(cusolver_handle_);
+    cublasDestroy(cublas_handle_);
+    return false;
+  }
+  if (cusolverDnSetStream(cusolver_handle_, stream_) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnSetStream failed.";
+    cusolverDnDestroy(cusolver_handle_);
+    cudaStreamDestroy(stream_);
+    cublasDestroy(cublas_handle_);
+    return false;
+  }
+  if (cublasSetStream(cublas_handle_, stream_) != CUBLAS_STATUS_SUCCESS) {
+    *message = "cuBLAS::cublasSetStream failed.";
+    cusolverDnDestroy(cusolver_handle_);
+    cublasDestroy(cublas_handle_);
+    cudaStreamDestroy(stream_);
+    return false;
+  }
+  error_.Reserve(1);
+  *message = "CUDADenseQR::Init Success.";
+  return true;
+}
+
+CUDADenseQR::~CUDADenseQR() {
+  if (cublas_handle_ != nullptr) {
+    CHECK_EQ(cublasDestroy(cublas_handle_), CUBLAS_STATUS_SUCCESS);
+  }
+  if (cusolver_handle_ != nullptr) {
+    CHECK_EQ(cusolverDnDestroy(cusolver_handle_), CUSOLVER_STATUS_SUCCESS);
+  }
+  if (stream_ != nullptr) {
+    CHECK_EQ(cudaStreamDestroy(stream_), cudaSuccess);
+  }
+}
+
+LinearSolverTerminationType CUDADenseQR::Factorize(
+    int num_rows, int num_cols, double* lhs, std::string* message) {
+  factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  lhs_.Reserve(num_rows * num_cols);
+  tau_.Reserve(std::min(num_rows, num_cols));
+  num_rows_ = num_rows;
+  num_cols_ = num_cols;
+  lhs_.CopyToGpuAsync(lhs, num_rows * num_cols, stream_);
+  int device_workspace_size = 0;
+  if (cusolverDnDgeqrf_bufferSize(cusolver_handle_,
+                                  num_rows,
+                                  num_cols,
+                                  lhs_.data(),
+                                  num_rows,
+                                  &device_workspace_size) !=
+      CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDgeqrf_bufferSize failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  device_workspace_.Reserve(device_workspace_size);
+  if (cusolverDnDgeqrf(cusolver_handle_,
+                       num_rows,
+                       num_cols,
+                       lhs_.data(),
+                       num_rows,
+                       tau_.data(),
+                       reinterpret_cast<double*>(device_workspace_.data()),
+                       device_workspace_.size(),
+                       error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDgeqrf 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;
+  error_.CopyToHost(&error, 1);
+  if (error < 0) {
+    LOG(FATAL) << "Congratulations, you found a bug in Ceres - "
+               << "please report it. "
+               << "cuSolverDN::cusolverDnDgeqrf 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;
+  }
+
+  *message = "Success";
+  factorize_result_ = LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+  return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+}
+
+LinearSolverTerminationType CUDADenseQR::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_;
+  }
+  rhs_.CopyToGpuAsync(rhs, num_rows_, stream_);
+  int device_workspace_size = 0;
+  if (cusolverDnDormqr_bufferSize(cusolver_handle_,
+                                  CUBLAS_SIDE_LEFT,
+                                  CUBLAS_OP_T,
+                                  num_rows_,
+                                  1,
+                                  num_cols_,
+                                  lhs_.data(),
+                                  num_rows_,
+                                  tau_.data(),
+                                  rhs_.data(),
+                                  num_rows_,
+                                  &device_workspace_size) !=
+      CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDormqr_bufferSize failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  device_workspace_.Reserve(device_workspace_size);
+  // Compute rhs = Q^T * rhs, assuming that lhs has already been factorized.
+  // The result of factorization would have stored Q in a packed form in lhs_.
+  if (cusolverDnDormqr(cusolver_handle_,
+                       CUBLAS_SIDE_LEFT,
+                       CUBLAS_OP_T,
+                       num_rows_,
+                       1,
+                       num_cols_,
+                       lhs_.data(),
+                       num_rows_,
+                       tau_.data(),
+                       rhs_.data(),
+                       num_rows_,
+                       reinterpret_cast<double*>(device_workspace_.data()),
+                       device_workspace_.size(),
+                       error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDormqr failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  int error = 0;
+  error_.CopyToHost(&error, 1);
+  if (error < 0) {
+    LOG(FATAL) << "Congratulations, you found a bug in Ceres. "
+               << "Please report it."
+               << "cuSolverDN::cusolverDnDormqr fatal error. "
+               << "Argument: " << -error << " is invalid.";
+  }
+  // Compute the solution vector as x = R \ (Q^T * rhs). Since the previous step
+  // replaced rhs by (Q^T * rhs), this is just x = R \ rhs.
+  if (cublasDtrsv(cublas_handle_,
+                  CUBLAS_FILL_MODE_UPPER,
+                  CUBLAS_OP_N,
+                  CUBLAS_DIAG_NON_UNIT,
+                  num_cols_,
+                  lhs_.data(),
+                  num_rows_,
+                  rhs_.data(),
+                  1) != CUBLAS_STATUS_SUCCESS) {
+    *message = "cuBLAS::cublasDtrsv failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_) != cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::LINEAR_SOLVER_FATAL_ERROR;
+  }
+  rhs_.CopyToHost(solution, num_cols_);
+  *message = "Success";
+  return LinearSolverTerminationType::LINEAR_SOLVER_SUCCESS;
+}
+
+std::unique_ptr<CUDADenseQR> CUDADenseQR::Create(
+    const LinearSolver::Options& options) {
+  if (options.dense_linear_algebra_library_type != CUDA) {
+    // The user called the wrong factory method.
+    return nullptr;
+  }
+  auto cuda_dense_qr =
+      std::unique_ptr<CUDADenseQR>(new CUDADenseQR());
+  std::string cuda_error;
+  if (cuda_dense_qr->Init(&cuda_error)) {
+    return cuda_dense_qr;
+  }
+  // Initialization failed, destroy the object (done automatically) and return a
+  // nullptr.
+  LOG(ERROR) << "CUDADenseQR::Init failed: " << cuda_error;
+  return nullptr;
+}
+
+CUDADenseQR::CUDADenseQR()  = default;
+
+#endif  // CERES_NO_CUDA
+
+
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/dense_qr.h b/internal/ceres/dense_qr.h
index aafb067..8bcccf0 100644
--- a/internal/ceres/dense_qr.h
+++ b/internal/ceres/dense_qr.h
@@ -45,6 +45,12 @@
 #include "ceres/internal/export.h"
 #include "ceres/linear_solver.h"
 #include "glog/logging.h"
+#ifndef CERES_NO_CUDA
+#include "ceres/cuda_buffer.h"
+#include "cuda_runtime.h"
+#include "cublas_v2.h"
+#include "cusolverDn.h"
+#endif  // CERES_NO_CUDA
 
 namespace ceres {
 namespace internal {
@@ -137,6 +143,65 @@
 };
 #endif  // CERES_NO_LAPACK
 
+#ifndef CERES_NO_CUDA
+// Implementation of DenseQR using the 32-bit cuSolverDn interface. A
+// requirement for using this solver is that the lhs must not be rank deficient.
+// This is because cuSolverDn does not implement the singularity-checking
+// wrapper trtrs, hence this solver directly uses trsv from CUBLAS for the
+// backsubstitution.
+class CERES_NO_EXPORT CUDADenseQR : public DenseQR {
+ public:
+  static std::unique_ptr<CUDADenseQR> Create(
+      const LinearSolver::Options& options);
+  ~CUDADenseQR() override;
+  CUDADenseQR(const CUDADenseQR&) = delete;
+  CUDADenseQR& operator=(const CUDADenseQR&) = delete;
+  LinearSolverTerminationType Factorize(int num_rows,
+                                        int num_cols,
+                                        double* lhs,
+                                        std::string* message) override;
+  LinearSolverTerminationType Solve(const double* rhs,
+                                    double* solution,
+                                    std::string* message) override;
+
+ private:
+  CUDADenseQR();
+  // 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;
+  // Handle to cuBLAS context.
+  cublasHandle_t cublas_handle_ = nullptr;
+  // CUDA device stream.
+  cudaStream_t stream_ = nullptr;
+  // Number of rowns in the A matrix, to be cached between calls to *Factorize
+  // and *Solve.
+  size_t num_rows_ = 0;
+  // 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_;
+  // GPU memory allocated for the TAU matrix (scaling of householder vectors).
+  CudaBuffer<double> tau_;
+  // 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;
+};
+
+#endif  // CERES_NO_CUDA
+
 }  // namespace internal
 }  // namespace ceres
 
diff --git a/internal/ceres/dense_qr_test.cc b/internal/ceres/dense_qr_test.cc
index ac1258d..d9d307e 100644
--- a/internal/ceres/dense_qr_test.cc
+++ b/internal/ceres/dense_qr_test.cc
@@ -33,6 +33,7 @@
 #include <memory>
 #include <numeric>
 #include <string>
+#include <tuple>
 #include <vector>
 
 #include "Eigen/Dense"
@@ -67,6 +68,7 @@
 
   LinearSolver::Options options;
   options.dense_linear_algebra_library_type = GetParam();
+  const double kEpsilon = std::numeric_limits<double>::epsilon() * 1.5e4;
   std::unique_ptr<DenseQR> dense_qr = DenseQR::Create(options);
 
   const int kNumTrials = 10;
@@ -83,7 +85,6 @@
         Vector x = VectorType::Random(num_cols);
         Vector rhs = lhs * x;
         Vector actual = Vector::Random(num_cols);
-
         LinearSolver::Summary summary;
         summary.termination_type = dense_qr->FactorAndSolve(num_rows,
                                                             num_cols,
@@ -91,10 +92,8 @@
                                                             rhs.data(),
                                                             actual.data(),
                                                             &summary.message);
-        EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
-        EXPECT_NEAR((x - actual).norm() / x.norm(),
-                    0.0,
-                    std::numeric_limits<double>::epsilon() * 400)
+        ASSERT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
+        ASSERT_NEAR((x - actual).norm() / x.norm(), 0.0, kEpsilon)
             << "\nexpected: " << x.transpose()
             << "\nactual  : " << actual.transpose();
       }
@@ -102,17 +101,26 @@
   }
 }
 
+namespace {
+
+// NOTE: preprocessor directives in a macro are not standard conforming
+decltype(auto) MakeValues() {
+  return ::testing::Values(
+      EIGEN
 #ifndef CERES_NO_LAPACK
-INSTANTIATE_TEST_SUITE_P(_,
-                         DenseQRTest,
-                         ::testing::Values(EIGEN, LAPACK),
-                         ParamInfoToString);
-#else
-INSTANTIATE_TEST_SUITE_P(_,
-                         DenseQRTest,
-                         ::testing::Values(EIGEN),
-                         ParamInfoToString);
+      ,
+      LAPACK
 #endif
+#ifndef CERES_NO_CUDA
+      ,
+      CUDA
+#endif
+  );
+}
+
+}  // namespace
+
+INSTANTIATE_TEST_SUITE_P(_, DenseQRTest, MakeValues(), ParamInfoToString);
 
 }  // namespace internal
 }  // namespace ceres