diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 27ea80e..b05c3a2 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -81,8 +81,9 @@
               "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 @@
             "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 @@
   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 1b07ab3..0a110e2 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -137,7 +137,7 @@
   // 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 @@
   // 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 99d5375..f4f9e0f 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() = 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 @@
 
   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 @@
   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<CudaVector> {
+ 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<CudaVector> {
+ 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> 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<CudaCgnrSolver>(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<CudaSparseMatrix>(options_.context, A);
+    b_ = std::make_unique<CudaVector>(options_.context, A.num_rows());
+    x_ = std::make_unique<CudaVector>(options_.context, A.num_cols());
+    Atb_ = std::make_unique<CudaVector>(options_.context, A.num_cols());
+    Ax_ = std::make_unique<CudaVector>(options_.context, A.num_rows());
+    D_ = std::make_unique<CudaVector>(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 6982296..a971c81 100644
--- a/internal/ceres/cgnr_solver.h
+++ b/internal/ceres/cgnr_solver.h
@@ -33,6 +33,8 @@
 
 #include <memory>
 
+#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 @@
   const LinearSolver::Options options_;
   std::unique_ptr<Preconditioner> 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<CudaCgnrSolver> 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<CudaSparseMatrix> A_;
+  std::unique_ptr<CudaVector> b_;
+  std::unique_ptr<CudaVector> x_;
+  std::unique_ptr<CudaVector> Atb_;
+  std::unique_ptr<CudaVector> Ax_;
+  std::unique_ptr<CudaVector> 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 f2eea87..0c3581b 100644
--- a/internal/ceres/conjugate_gradients_solver.h
+++ b/internal/ceres/conjugate_gradients_solver.h
@@ -100,25 +100,25 @@
 // 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 <typename DenseVectorType>
 LinearSolver::Summary ConjugateGradientsSolver(
     const ConjugateGradientsSolverOptions options,
     ConjugateGradientsLinearOperator<DenseVectorType>& lhs,
     const DenseVectorType& rhs,
     ConjugateGradientsLinearOperator<DenseVectorType>& 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 @@
     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 @@
     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 a01dfc9..e59eba7 100644
--- a/internal/ceres/conjugate_gradients_solver_test.cc
+++ b/internal/ceres/conjugate_gradients_solver_test.cc
@@ -74,9 +74,10 @@
   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 @@
   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 c7b11ce..1483923 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 @@
   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 @@
   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 @@
   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 @@
   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 @@
   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 @@
   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<DenseCholesky> dense_cholesky =
       CUDADenseCholesky::Create(options);
@@ -175,12 +186,20 @@
   {
     // 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 @@
   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 @@
   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 @@
   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 @@
 
 #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 9eb5d4d..cc0dbd6 100644
--- a/internal/ceres/cuda_dense_qr_test.cc
+++ b/internal/ceres/cuda_dense_qr_test.cc
@@ -41,6 +41,10 @@
 
 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 @@
   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 @@
   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 @@
   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 @@
   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<DenseQR> 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 244de74..af42e73 100644
--- a/internal/ceres/dense_cholesky.cc
+++ b/internal/ceres/dense_cholesky.cc
@@ -348,9 +348,8 @@
 #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 @@
 
 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 1c0df56..34ff4f3 100644
--- a/internal/ceres/dense_cholesky_test.cc
+++ b/internal/ceres/dense_cholesky_test.cc
@@ -76,6 +76,8 @@
   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 acef750..f562d97 100644
--- a/internal/ceres/dense_qr_test.cc
+++ b/internal/ceres/dense_qr_test.cc
@@ -68,6 +68,8 @@
   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<double>::epsilon() * 1.5e4;
   std::unique_ptr<DenseQR> dense_qr = DenseQR::Create(options);
diff --git a/internal/ceres/evaluator.cc b/internal/ceres/evaluator.cc
index f44796a..91cf221 100644
--- a/internal/ceres/evaluator.cc
+++ b/internal/ceres/evaluator.cc
@@ -64,10 +64,17 @@
     case DENSE_SCHUR:
     case SPARSE_SCHUR:
     case ITERATIVE_SCHUR:
-    case CGNR:
-      return std::make_unique<
-          ProgramEvaluator<BlockEvaluatePreparer, BlockJacobianWriter>>(
-          options, program);
+    case CGNR: {
+      if (options.sparse_linear_algebra_library_type == CUDA_SPARSE) {
+        return std::make_unique<
+          ProgramEvaluator<ScratchEvaluatePreparer,
+                           CompressedRowJacobianWriter>>(options, program);
+      } else {
+        return std::make_unique<
+            ProgramEvaluator<BlockEvaluatePreparer, BlockJacobianWriter>>(
+            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 68a4fb2..7ea523c 100644
--- a/internal/ceres/evaluator.h
+++ b/internal/ceres/evaluator.h
@@ -65,6 +65,7 @@
     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 607502e..0463afb 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -129,6 +129,9 @@
   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 @@
                                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 d9ef150..2e3fe4b 100644
--- a/internal/ceres/linear_solver.cc
+++ b/internal/ceres/linear_solver.cc
@@ -76,8 +76,15 @@
   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<CgnrSolver>(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 ba96e1c..631862c 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -217,7 +217,14 @@
   }
 }
 
-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 @@
   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<Vector>(
       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 4fde69c..cac96b5 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/internal/ceres/schur_complement_solver.h
@@ -188,7 +188,7 @@
   std::unique_ptr<SparseCholesky> sparse_cholesky_;
   std::unique_ptr<BlockRandomAccessDiagonalMatrix> 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 b772f3f..0e0790e 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -327,7 +327,8 @@
     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 @@
     }
   }
 
+  // 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 @@
   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 @@
   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 @@
       }
     }
 
-    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 b4b2d34..92c1ddc 100644
--- a/internal/ceres/solver_test.cc
+++ b/internal/ceres/solver_test.cc
@@ -878,6 +878,20 @@
   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 @@
   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 @@
     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 d27a472..84e783f 100644
--- a/internal/ceres/trust_region_preprocessor.cc
+++ b/internal/ceres/trust_region_preprocessor.cc
@@ -261,6 +261,8 @@
   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 1d514de..3de91ed 100644
--- a/internal/ceres/types.cc
+++ b/internal/ceres/types.cc
@@ -108,6 +108,7 @@
     CASESTR(SUITE_SPARSE);
     CASESTR(EIGEN_SPARSE);
     CASESTR(ACCELERATE_SPARSE);
+    CASESTR(CUDA_SPARSE);
     CASESTR(NO_SPARSE);
     default:
       return "UNKNOWN";
@@ -120,6 +121,7 @@
   STRENUM(SUITE_SPARSE);
   STRENUM(EIGEN_SPARSE);
   STRENUM(ACCELERATE_SPARSE);
+  STRENUM(CUDA_SPARSE);
   STRENUM(NO_SPARSE);
   return false;
 }
@@ -420,6 +422,14 @@
 #endif
   }
 
+  if (type == CUDA_SPARSE) {
+#ifdef CERES_NO_CUDA
+    return false;
+#else
+    return true;
+#endif
+  }
+
   if (type == NO_SPARSE) {
     return true;
   }