diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8d6094c..82014a9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -145,7 +145,9 @@
   option(ENABLE_BITCODE
     "Enable bitcode for iOS builds (disables inline optimizations for Eigen)." OFF)
 endif()
-option(CUDA "Enable use of CUDA linear algebra solvers." ON)
+# We can't have an option called 'CUDA' since that is a reserved word -- a
+# language definition.
+option(USE_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
@@ -228,23 +230,23 @@
   endif (EIGENSPARSE)
 endif (Eigen3_FOUND)
 
-if (CUDA)
+if (USE_CUDA)
   find_package(CUDA QUIET)
   if (CUDA_FOUND)
     message("-- Found CUDA version ${CUDA_VERSION}: "
         "${CUDA_LIBRARIES};"
         "${CUDA_cusolver_LIBRARY};"
         "${CUDA_cusparse_LIBRARY}")
+    enable_language(CUDA)
   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)
+else (USE_CUDA)
   message("-- Building without CUDA.")
   list(APPEND CERES_COMPILE_OPTIONS CERES_NO_CUDA)
-endif (CUDA)
-
+endif (USE_CUDA)
 if (LAPACK)
   find_package(LAPACK QUIET)
   if (LAPACK_FOUND)
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 5e6b0a7..d6bae75 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -118,11 +118,12 @@
 
 - `CUDA <https://developer.nvidia.com/cuda-toolkit>`_ If you have an
   NVIDIA GPU then Ceres Solver can use it accelerate the solution of
-  the Gauss-Newton linear systems using ``CUDA``. Currently this
-  support is limited to using the dense linear solvers that ship with
-  ``CUDA``. As a result GPU acceleration can be used to speed up
+  the Gauss-Newton linear systems using the CMake flag ``USE_CUDA``.
+  Currently this support is limited to using the dense linear solvers that ship
+  with ``CUDA``. As a result GPU acceleration can be used to speed up
   ``DENSE_QR``, ``DENSE_NORMAL_CHOLESKY`` and
-  ``DENSE_SCHUR``. **Optional**.
+  ``DENSE_SCHUR``. This also enables ``CUDA`` mixed precision solves
+  for ``DENSE_NORMAL_CHOLESKY`` and ``DENSE_SCHUR``.  **Optional**.
 
 .. _section-linux:
 
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 097276b..3b15318 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -573,6 +573,8 @@
 
     // TODO(sameeragarwal): Further expand the documentation for the
     // following two options.
+    // TODO(joydeepbiswas): Update the documentation for the mixed precision
+    // option with CUDA.
 
     // NOTE1: EXPERIMENTAL FEATURE, UNDER DEVELOPMENT, USE AT YOUR OWN RISK.
     //
@@ -996,6 +998,8 @@
         SPARSE_NORMAL_CHOLESKY;
 #endif
 
+    bool mixed_precision_solves_used = false;
+
     LinearSolverOrderingType linear_solver_ordering_type;
 
     // Size of the elimination groups given by the user as hints to
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 4e5008b..ae22a9d 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -157,6 +157,10 @@
         COMMAND ${CUDA_TOOLKIT_ROOT_DIR}/bin/cuda-memcheck --leak-check full
             $<TARGET_FILE:cuda_dense_cholesky_test>)
   endif (BUILD_TESTING)
+  set_source_files_properties(ceres_cuda_kernels.cu PROPERTIES LANGUAGE CUDA)
+  add_library(ceres_cuda_kernels ceres_cuda_kernels.cu)
+  target_compile_features(ceres_cuda_kernels PRIVATE cxx_std_14)
+  list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ceres_cuda_kernels)
 endif (CUDA_FOUND)
 
 if (LAPACK_FOUND)
@@ -387,7 +391,7 @@
   list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES_INCLUDE_DIRS
     ${AccelerateSparse_INCLUDE_DIRS})
 endif()
-if (CUDA)
+if (USE_CUDA)
   list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES_INCLUDE_DIRS
     ${CUDA_INCLUDE_DIRS})
 endif()
@@ -403,6 +407,12 @@
 generate_export_header(ceres EXPORT_FILE_NAME
   ${Ceres_BINARY_DIR}/${CMAKE_INSTALL_INCLUDEDIR}/ceres/internal/export.h)
 
+install(TARGETS ceres_cuda_kernels
+        EXPORT  CeresExport
+        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
+
 install(TARGETS ceres
         EXPORT  CeresExport
         RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
diff --git a/internal/ceres/ceres_cuda_kernels.cu b/internal/ceres/ceres_cuda_kernels.cu
new file mode 100644
index 0000000..2fd76b5
--- /dev/null
+++ b/internal/ceres/ceres_cuda_kernels.cu
@@ -0,0 +1,109 @@
+// 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 "cuda_runtime.h"
+
+namespace ceres::internal {
+
+// As the CUDA Toolkit documentation says, "although arbitrary in this case, is
+// a common choice". This is determined by the warp size, max block size, and
+// multiprocessor sizes of recent GPUs. For complex kernels with significant
+// register usage and unusual memory patterns, the occupancy calculator API
+// might provide better performance. See "Occupancy Calculator" under the CUDA
+// toolkit documentation.
+constexpr int kCudaBlockSize = 256;
+
+template<typename SrcType, typename DstType>
+__global__ void TypeConversionKernel(const SrcType* __restrict__ input,
+                                     DstType* __restrict__ output,
+                                     const int size) {
+  const int i = blockIdx.x * blockDim.x + threadIdx.x;
+  if (i < size) {
+    output[i] = static_cast<DstType>(input[i]);
+  }
+}
+
+void CudaFP64ToFP32(const double* input,
+                    float* output,
+                    const int size,
+                    cudaStream_t stream) {
+  const int num_blocks = (size + kCudaBlockSize - 1) / kCudaBlockSize;
+  TypeConversionKernel<double, float>
+      <<<num_blocks, kCudaBlockSize, 0, stream>>>(input, output, size);
+}
+
+void CudaFP32ToFP64(const float* input,
+                    double* output,
+                    const int size,
+                    cudaStream_t stream) {
+  const int num_blocks = (size + kCudaBlockSize - 1) / kCudaBlockSize;
+  TypeConversionKernel<float, double>
+      <<<num_blocks, kCudaBlockSize, 0, stream>>>(input, output, size);
+}
+
+template<typename T>
+__global__ void SetZeroKernel(T* __restrict__ output, const int size) {
+  const int i = blockIdx.x * blockDim.x + threadIdx.x;
+  if (i < size) {
+    output[i] = T(0.0);
+  }
+}
+
+void CudaSetZeroFP32(float* output, const int size, cudaStream_t stream) {
+  const int num_blocks = (size + kCudaBlockSize - 1) / kCudaBlockSize;
+  SetZeroKernel<float><<<num_blocks, kCudaBlockSize, 0, stream>>>(output, size);
+}
+
+void CudaSetZeroFP64(double* output, const int size, cudaStream_t stream) {
+  const int num_blocks = (size + kCudaBlockSize - 1) / kCudaBlockSize;
+  SetZeroKernel<double><<<num_blocks, kCudaBlockSize, 0, stream>>>(
+      output, size);
+}
+
+template <typename SrcType, typename DstType>
+__global__ void XPlusEqualsYKernel(DstType* __restrict__ x,
+                                   const SrcType* __restrict__ y,
+                                  const int size) {
+  const int i = blockIdx.x * blockDim.x + threadIdx.x;
+  if (i < size) {
+    x[i] = x[i] + DstType(y[i]);
+  }
+}
+
+void CudaDsxpy(double* x,
+                float* y,
+                const int size,
+                cudaStream_t stream) {
+  const int num_blocks = (size + kCudaBlockSize - 1) / kCudaBlockSize;
+  XPlusEqualsYKernel<float, double>
+      <<<num_blocks, kCudaBlockSize, 0, stream>>>(x, y, size);
+}
+
+} // namespace ceres_cuda_kernels
\ No newline at end of file
diff --git a/internal/ceres/ceres_cuda_kernels.h b/internal/ceres/ceres_cuda_kernels.h
new file mode 100644
index 0000000..989399f
--- /dev/null
+++ b/internal/ceres/ceres_cuda_kernels.h
@@ -0,0 +1,66 @@
+// 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 "ceres/internal/config.h"
+
+#ifndef CERES_NO_CUDA
+
+#include "cuda_runtime.h"
+
+namespace ceres::internal {
+
+// Convert an array of double (FP64) values to float (FP32). Both arrays must
+// already be on GPU memory.
+void CudaFP64ToFP32(const double* input,
+                    float* output,
+                    const int size,
+                    cudaStream_t stream);
+
+// Convert an array of float (FP32) values to double (FP64). Both arrays must
+// already be on GPU memory.
+void CudaFP32ToFP64(const float* input,
+                    double* output,
+                    const int size,
+                    cudaStream_t stream);
+
+// Set all elements of the array to the FP32 value 0. The array must be in GPU
+// memory.
+void CudaSetZeroFP32(float* output, const int size, cudaStream_t stream);
+
+// Set all elements of the array to the FP64 value 0. The array must be in GPU
+// memory.
+void CudaSetZeroFP64(double* output, const int size, cudaStream_t stream);
+
+// Compute x = x + double(y). Input array is float (FP32), output array is
+// double (FP64). Both arrays must already be on GPU memory.
+void CudaDsxpy(double* x, float* y, const int size, cudaStream_t stream);
+}  // namespace ceres_cuda_kernels
+
+#endif  // CERES_NO_CUDA
\ No newline at end of file
diff --git a/internal/ceres/cuda_buffer.h b/internal/ceres/cuda_buffer.h
index a1cf784..8868e1a 100644
--- a/internal/ceres/cuda_buffer.h
+++ b/internal/ceres/cuda_buffer.h
@@ -81,6 +81,14 @@
              cudaSuccess);
   }
 
+  // Perform an asynchronous copy from GPU memory using the stream provided.
+  void CopyFromGpuAsync(const T* data, const size_t size, cudaStream_t stream) {
+    Reserve(size);
+    CHECK_EQ(cudaMemcpyAsync(
+        data_, data, size * sizeof(T), cudaMemcpyDeviceToDevice, 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.
@@ -95,6 +103,7 @@
   }
 
   T* data() { return data_; }
+  const T* data() const { return data_; }
   size_t size() const { return size_; }
 
  private:
diff --git a/internal/ceres/cuda_dense_cholesky_test.cc b/internal/ceres/cuda_dense_cholesky_test.cc
index 13dc34b..b9acc99 100644
--- a/internal/ceres/cuda_dense_cholesky_test.cc
+++ b/internal/ceres/cuda_dense_cholesky_test.cc
@@ -165,8 +165,140 @@
     summary.termination_type = dense_cholesky->FactorAndSolve(
         kNumCols, lhs.data(), rhs.data(), x_computed.data(), &summary.message);
     ASSERT_EQ(summary.termination_type, LinearSolverTerminationType::SUCCESS);
+    static const double kEpsilon = std::numeric_limits<double>::epsilon() * 2e5;
     ASSERT_NEAR(
-        (x_computed - x_expected).norm() / x_expected.norm(), 0.0, 1e-10);
+        (x_computed - x_expected).norm() / x_expected.norm(), 0.0, kEpsilon);
+  }
+}
+
+TEST(CUDADenseCholeskyMixedPrecision, InvalidOptionsOnCreate) {
+  {
+    // Did not ask for CUDA, and did not ask for mixed precision.
+    LinearSolver::Options options;
+    auto solver = CUDADenseCholeskyMixedPrecision::Create(options);
+    ASSERT_EQ(solver, nullptr);
+  }
+  {
+    // Asked for CUDA, but did not ask for mixed precision.
+    LinearSolver::Options options;
+    options.dense_linear_algebra_library_type = ceres::CUDA;
+    auto solver = CUDADenseCholeskyMixedPrecision::Create(options);
+    ASSERT_EQ(solver, nullptr);
+  }
+}
+
+// Tests the CUDA Cholesky solver with a simple 4x4 matrix.
+TEST(CUDADenseCholeskyMixedPrecision, Cholesky4x4Matrix1Step) {
+  Eigen::Matrix4d A;
+  // clang-format off
+  // A common test Cholesky decomposition test matrix, see :
+  // https://en.wikipedia.org/w/index.php?title=Cholesky_decomposition&oldid=1080607368#Example
+  A <<  4,  12, -16, 0,
+       12,  37, -43, 0,
+      -16, -43,  98, 0,
+        0,   0,   0, 1;
+  // clang-format on
+
+  const Eigen::Vector4d b = Eigen::Vector4d::Ones();
+  LinearSolver::Options options;
+  options.max_num_refinement_iterations = 0;
+  ContextImpl context;
+  options.context = &context;
+  options.dense_linear_algebra_library_type = CUDA;
+  options.use_mixed_precision_solves = true;
+  auto solver = CUDADenseCholeskyMixedPrecision::Create(options);
+  ASSERT_NE(solver, nullptr);
+  std::string error_string;
+  ASSERT_EQ(solver->Factorize(A.cols(), A.data(), &error_string),
+            LinearSolverTerminationType::SUCCESS);
+  Eigen::Vector4d x = Eigen::Vector4d::Zero();
+  ASSERT_EQ(solver->Solve(b.data(), x.data(), &error_string),
+            LinearSolverTerminationType::SUCCESS);
+  // A single step of the mixed precision solver will be equivalent to solving
+  // in low precision (FP32). Hence the tolerance is defined w.r.t. FP32 epsilon
+  // instead of FP64 epsilon.
+  static const double kEpsilon = std::numeric_limits<float>::epsilon() * 15.0f;
+  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 Cholesky solver with a simple 4x4 matrix.
+TEST(CUDADenseCholeskyMixedPrecision, Cholesky4x4Matrix4Steps) {
+  Eigen::Matrix4d A;
+  // clang-format off
+  A <<  4,  12, -16, 0,
+       12,  37, -43, 0,
+      -16, -43,  98, 0,
+        0,   0,   0, 1;
+  // clang-format on
+
+  const Eigen::Vector4d b = Eigen::Vector4d::Ones();
+  LinearSolver::Options options;
+  options.max_num_refinement_iterations = 3;
+  ContextImpl context;
+  options.context = &context;
+  options.dense_linear_algebra_library_type = CUDA;
+  options.use_mixed_precision_solves = true;
+  auto solver = CUDADenseCholeskyMixedPrecision::Create(options);
+  ASSERT_NE(solver, nullptr);
+  std::string error_string;
+  ASSERT_EQ(solver->Factorize(A.cols(), A.data(), &error_string),
+            LinearSolverTerminationType::SUCCESS);
+  Eigen::Vector4d x = Eigen::Vector4d::Zero();
+  ASSERT_EQ(solver->Solve(b.data(), x.data(), &error_string),
+            LinearSolverTerminationType::SUCCESS);
+  // The error does not reduce beyond four iterations, and stagnates at this
+  // level of precision.
+  static const double kEpsilon = std::numeric_limits<double>::epsilon() * 3e2;
+  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);
+}
+
+
+
+TEST(CUDADenseCholeskyMixedPrecision, Randomized1600x1600Tests) {
+  const int kNumCols = 1600;
+  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;
+  ContextImpl context;
+  options.context = &context;
+  options.dense_linear_algebra_library_type = ceres::CUDA;
+  options.use_mixed_precision_solves = true;
+  options.max_num_refinement_iterations = 20;
+  std::unique_ptr<CUDADenseCholeskyMixedPrecision> dense_cholesky =
+      CUDADenseCholeskyMixedPrecision::Create(options);
+
+  const int kNumTrials = 20;
+  for (int i = 0; i < kNumTrials; ++i) {
+    LhsType lhs = LhsType::Random(kNumCols, kNumCols);
+    lhs = lhs.transpose() * lhs;
+    lhs += 1e-3 * LhsType::Identity(kNumCols, 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(), kNumCols);
+    EXPECT_EQ(lhs.cols(), kNumCols);
+    EXPECT_EQ(rhs.rows(), kNumCols);
+    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_cholesky->FactorAndSolve(
+        kNumCols, lhs.data(), rhs.data(), x_computed.data(), &summary.message);
+    ASSERT_EQ(summary.termination_type, LinearSolverTerminationType::SUCCESS);
+    static const double kEpsilon = std::numeric_limits<double>::epsilon() * 1e6;
+    ASSERT_NEAR(
+        (x_computed - x_expected).norm() / x_expected.norm(), 0.0, kEpsilon);
   }
 }
 
diff --git a/internal/ceres/dense_cholesky.cc b/internal/ceres/dense_cholesky.cc
index a10f311..e2c5036 100644
--- a/internal/ceres/dense_cholesky.cc
+++ b/internal/ceres/dense_cholesky.cc
@@ -38,6 +38,7 @@
 #include "ceres/internal/config.h"
 
 #ifndef CERES_NO_CUDA
+#include "ceres/ceres_cuda_kernels.h"
 #include "ceres/context_impl.h"
 #include "cuda_runtime.h"
 #include "cusolverDn.h"
@@ -68,12 +69,16 @@
   std::unique_ptr<DenseCholesky> dense_cholesky;
 
   switch (options.dense_linear_algebra_library_type) {
-    case EIGEN:
+    case EIGEN: {
+      // Eigen mixed precision solver not yet implemented.
+      if (options.use_mixed_precision_solves) return nullptr;
       dense_cholesky = std::make_unique<EigenDenseCholesky>();
-      break;
+    } break;
 
     case LAPACK:
 #ifndef CERES_NO_LAPACK
+      // LAPACK mixed precision solver not yet implemented.
+      if (options.use_mixed_precision_solves) return nullptr;
       dense_cholesky = std::make_unique<LAPACKDenseCholesky>();
       break;
 #else
@@ -82,7 +87,11 @@
 
     case CUDA:
 #ifndef CERES_NO_CUDA
-      dense_cholesky = CUDADenseCholesky::Create(options);
+      if (options.use_mixed_precision_solves) {
+        dense_cholesky = CUDADenseCholeskyMixedPrecision::Create(options);
+      } else {
+        dense_cholesky = CUDADenseCholesky::Create(options);
+      }
       break;
 #else
       LOG(FATAL) << "Ceres was compiled without support for CUDA.";
@@ -320,6 +329,210 @@
   return nullptr;
 }
 
+std::unique_ptr<CUDADenseCholeskyMixedPrecision>
+    CUDADenseCholeskyMixedPrecision::Create(
+    const LinearSolver::Options& options) {
+  if (options.dense_linear_algebra_library_type != CUDA ||
+      !options.use_mixed_precision_solves) {
+    // The user called the wrong factory method.
+    return nullptr;
+  }
+  auto solver = std::unique_ptr<CUDADenseCholeskyMixedPrecision>(
+      new CUDADenseCholeskyMixedPrecision());
+  std::string cuda_error;
+  if (solver->Init(options, &cuda_error)) {
+    return solver;
+  }
+  LOG(ERROR) << "CUDADenseCholeskyMixedPrecision::Init failed: " << cuda_error;
+  return nullptr;
+}
+
+bool CUDADenseCholeskyMixedPrecision::Init(
+    const LinearSolver::Options& options, std::string* message) {
+  if (!options.context->InitCUDA(message)) {
+    return false;
+  }
+  cusolver_handle_ = options.context->cusolver_handle_;
+  cublas_handle_ = options.context->cublas_handle_;
+  stream_ = options.context->stream_;
+  error_.Reserve(1);
+  max_num_refinement_iterations_ = options.max_num_refinement_iterations;
+  *message = "CUDADenseCholeskyMixedPrecision::Init Success.";
+  return true;
+}
+
+LinearSolverTerminationType
+CUDADenseCholeskyMixedPrecision::CudaCholeskyFactorize(std::string* message) {
+  int device_workspace_size = 0;
+  if (cusolverDnSpotrf_bufferSize(cusolver_handle_,
+                                  CUBLAS_FILL_MODE_LOWER,
+                                  num_cols_,
+                                  lhs_fp32_.data(),
+                                  num_cols_,
+                                  &device_workspace_size) !=
+      CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnSpotrf_bufferSize failed.";
+    return LinearSolverTerminationType::FATAL_ERROR;
+  }
+  device_workspace_.Reserve(device_workspace_size);
+  if (cusolverDnSpotrf(cusolver_handle_,
+                       CUBLAS_FILL_MODE_LOWER,
+                       num_cols_,
+                       lhs_fp32_.data(),
+                       num_cols_,
+                       device_workspace_.data(),
+                       device_workspace_.size(),
+                       error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnSpotrf failed.";
+    return LinearSolverTerminationType::FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_) != cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::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::cusolverDnSpotrf 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::FATAL_ERROR;
+  }
+  if (error > 0) {
+    *message = StringPrintf(
+        "cuSolverDN::cusolverDnSpotrf numerical failure. "
+        "The leading minor of order %d is not positive definite.",
+        error);
+    factorize_result_ = LinearSolverTerminationType::FAILURE;
+    return LinearSolverTerminationType::FAILURE;
+  }
+  *message = "Success";
+  return LinearSolverTerminationType::SUCCESS;
+}
+
+LinearSolverTerminationType CUDADenseCholeskyMixedPrecision::CudaCholeskySolve(
+    std::string* message) {
+  CHECK_EQ(cudaMemcpyAsync(correction_fp32_.data(),
+                           residual_fp32_.data(),
+                           num_cols_ * sizeof(float),
+                           cudaMemcpyDeviceToDevice,
+                           stream_),
+           cudaSuccess);
+  if (cusolverDnSpotrs(cusolver_handle_,
+                       CUBLAS_FILL_MODE_LOWER,
+                       num_cols_,
+                       1,
+                       lhs_fp32_.data(),
+                       num_cols_,
+                       correction_fp32_.data(),
+                       num_cols_,
+                       error_.data()) != CUSOLVER_STATUS_SUCCESS) {
+    *message = "cuSolverDN::cusolverDnDpotrs failed.";
+    return LinearSolverTerminationType::FATAL_ERROR;
+  }
+  if (cudaDeviceSynchronize() != cudaSuccess ||
+      cudaStreamSynchronize(stream_) != cudaSuccess) {
+    *message = "Cuda device synchronization failed.";
+    return LinearSolverTerminationType::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::cusolverDnDpotrs fatal error. "
+               << "Argument: " << -error << " is invalid.";
+  }
+  *message = "Success";
+  return LinearSolverTerminationType::SUCCESS;
+}
+
+LinearSolverTerminationType CUDADenseCholeskyMixedPrecision::Factorize(
+    int num_cols,
+    double* lhs,
+    std::string* message) {
+  num_cols_ = num_cols;
+
+  // Copy fp64 version of lhs to GPU.
+  lhs_fp64_.Reserve(num_cols * num_cols);
+  lhs_fp64_.CopyToGpuAsync(lhs, num_cols * num_cols, stream_);
+
+  // Create an fp32 copy of lhs, lhs_fp32.
+  lhs_fp32_.Reserve(num_cols * num_cols);
+  CudaFP64ToFP32(
+      lhs_fp64_.data(), lhs_fp32_.data(), num_cols * num_cols, stream_);
+
+  // Factorize lhs_fp32.
+  factorize_result_ = CudaCholeskyFactorize(message);
+  return factorize_result_;
+}
+
+LinearSolverTerminationType CUDADenseCholeskyMixedPrecision::Solve(
+    const double* rhs,
+    double* solution,
+    std::string* message) {
+  // If factorization failed, return failure.
+  if (factorize_result_ != LinearSolverTerminationType::SUCCESS) {
+    *message = "Factorize did not complete successfully previously.";
+    return factorize_result_;
+  }
+
+  // Reserve memory for all arrays.
+  rhs_fp64_.Reserve(num_cols_);
+  x_fp64_.Reserve(num_cols_);
+  correction_fp32_.Reserve(num_cols_);
+  residual_fp32_.Reserve(num_cols_);
+  residual_fp64_.Reserve(num_cols_);
+
+  // Initialize x = 0.
+  CudaSetZeroFP64(x_fp64_.data(), num_cols_, stream_);
+
+  // Initialize residual = rhs.
+  rhs_fp64_.CopyToGpuAsync(rhs, num_cols_, stream_);
+  residual_fp64_.CopyFromGpuAsync(rhs_fp64_.data(), num_cols_, stream_);
+
+  for (int i = 0; i <= max_num_refinement_iterations_; ++i) {
+    // Cast residual from fp64 to fp32.
+    CudaFP64ToFP32(
+        residual_fp64_.data(), residual_fp32_.data(), num_cols_, stream_);
+    // [fp32] c = lhs^-1 * residual.
+    auto result = CudaCholeskySolve(message);
+    if (result != LinearSolverTerminationType::SUCCESS) {
+      return result;
+    }
+    // [fp64] x += c.
+    CudaDsxpy(
+        x_fp64_.data(), correction_fp32_.data(), num_cols_, stream_);
+    if (i < max_num_refinement_iterations_) {
+      // [fp64] residual = rhs - lhs * x
+      // This is done in two steps:
+      // 1. [fp64] residual = rhs
+      residual_fp64_.CopyFromGpuAsync(rhs_fp64_.data(), num_cols_, stream_);
+      // 2. [fp64] residual = residual - lhs * x
+      double alpha = -1.0;
+      double beta = 1.0;
+      cublasDsymv(cublas_handle_,
+                  CUBLAS_FILL_MODE_LOWER,
+                  num_cols_,
+                  &alpha,
+                  lhs_fp64_.data(),
+                  num_cols_,
+                  x_fp64_.data(),
+                  1,
+                  &beta,
+                  residual_fp64_.data(),
+                  1);
+    }
+  }
+  x_fp64_.CopyToHost(solution, num_cols_);
+  *message = "Success.";
+  return LinearSolverTerminationType::SUCCESS;
+}
+
 #endif  // CERES_NO_CUDA
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/dense_cholesky.h b/internal/ceres/dense_cholesky.h
index cc8642c..0593875 100644
--- a/internal/ceres/dense_cholesky.h
+++ b/internal/ceres/dense_cholesky.h
@@ -176,6 +176,78 @@
       LinearSolverTerminationType::FATAL_ERROR;
 };
 
+// A mixed-precision iterative refinement dense Cholesky solver using FP32 CUDA
+// Dense Cholesky for inner iterations, and FP64 outer refinements.
+// This class implements a modified version of the  "Classical iterative
+// refinement" (Algorithm 4.1) from the following paper:
+// Haidar, Azzam, Harun Bayraktar, Stanimire Tomov, Jack Dongarra, and Nicholas
+// J. Higham. "Mixed-precision iterative refinement using tensor cores on GPUs
+// to accelerate solution of linear systems." Proceedings of the Royal Society A
+// 476, no. 2243 (2020): 20200110.
+//
+// The three key modifications from Algorithm 4.1 in the paper are:
+// 1. We use Cholesky factorization instead of LU factorization since our A is
+//    symmetric positive definite.
+// 2. During the solution update, the up-cast and accumulation is performed in
+//    one step with a custom kernel.
+class CERES_NO_EXPORT CUDADenseCholeskyMixedPrecision final :
+    public DenseCholesky {
+ public:
+  static std::unique_ptr<CUDADenseCholeskyMixedPrecision> Create(
+      const LinearSolver::Options& options);
+  CUDADenseCholeskyMixedPrecision(
+      const CUDADenseCholeskyMixedPrecision&) = delete;
+  CUDADenseCholeskyMixedPrecision& operator=(
+      const CUDADenseCholeskyMixedPrecision&) = delete;
+  LinearSolverTerminationType Factorize(int num_cols,
+                                        double* lhs,
+                                        std::string* message) override;
+  LinearSolverTerminationType Solve(const double* rhs,
+                                    double* solution,
+                                    std::string* message) override;
+
+ private:
+  CUDADenseCholeskyMixedPrecision() = default;
+  // Helper function to wrap Cuda boilerplate needed to call Spotrf.
+  LinearSolverTerminationType CudaCholeskyFactorize(std::string* message);
+  // Helper function to wrap Cuda boilerplate needed to call Spotrs.
+  LinearSolverTerminationType CudaCholeskySolve(std::string* message);
+  // Picks up the cuSolverDN and cuStream handles from the context in the
+  // options, and the number of refinement iterations from the options. If
+  // the context is unable to initialize CUDA, returns false with a
+  // human-readable message indicating the reason.
+  bool Init(const LinearSolver::Options& options, std::string* message);
+
+  cusolverDnHandle_t cusolver_handle_ = nullptr;
+  cublasHandle_t cublas_handle_ = nullptr;
+  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;
+  CudaBuffer<double> lhs_fp64_;
+  CudaBuffer<double> rhs_fp64_;
+  CudaBuffer<float> lhs_fp32_;
+  // Scratch space for cuSOLVER on the GPU.
+  CudaBuffer<float> device_workspace_;
+  // Required for error handling with cuSOLVER.
+  CudaBuffer<int> error_;
+
+  // Solution to lhs * x = rhs.
+  CudaBuffer<double> x_fp64_;
+  // Incremental correction to x.
+  CudaBuffer<float> correction_fp32_;
+  // Residual to iterative refinement.
+  CudaBuffer<float> residual_fp32_;
+  CudaBuffer<double> residual_fp64_;
+
+  // Number of inner refinement iterations to perform.
+  int max_num_refinement_iterations_ = 0;
+  // Cache the result of Factorize to ensure that when Solve is called, the
+  // factorization of lhs is valid.
+  LinearSolverTerminationType factorize_result_ =
+      LinearSolverTerminationType::FATAL_ERROR;
+};
+
 #endif  // CERES_NO_CUDA
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/dense_cholesky_test.cc b/internal/ceres/dense_cholesky_test.cc
index 7866d5c..a425343 100644
--- a/internal/ceres/dense_cholesky_test.cc
+++ b/internal/ceres/dense_cholesky_test.cc
@@ -45,12 +45,20 @@
 
 namespace ceres::internal {
 
-using Param = DenseLinearAlgebraLibraryType;
+
+using Param =
+    ::testing::tuple<DenseLinearAlgebraLibraryType, bool>;
+constexpr bool kMixedPrecision = true;
+constexpr bool kFullPrecision = false;
 
 namespace {
 
 std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
-  return DenseLinearAlgebraLibraryTypeToString(info.param);
+  Param param = info.param;
+  std::stringstream ss;
+  ss << DenseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
+     << (::testing::get<1>(param) ? "MixedPrecision" : "FullPrecision");
+  return ss.str();
 }
 }  // namespace
 
@@ -67,14 +75,18 @@
   LinearSolver::Options options;
   ContextImpl context;
   options.context = &context;
-  options.dense_linear_algebra_library_type = GetParam();
+  options.dense_linear_algebra_library_type = ::testing::get<0>(GetParam());
+  options.use_mixed_precision_solves = ::testing::get<1>(GetParam());
+  const int kNumRefinementSteps = 4;
+  if (options.use_mixed_precision_solves) {
+    options.max_num_refinement_iterations = kNumRefinementSteps;
+  }
   std::unique_ptr<DenseCholesky> dense_cholesky =
       DenseCholesky::Create(options);
 
   const int kNumTrials = 10;
   const int kMinNumCols = 1;
   const int kMaxNumCols = 10;
-
   for (int num_cols = kMinNumCols; num_cols < kMaxNumCols; ++num_cols) {
     for (int trial = 0; trial < kNumTrials; ++trial) {
       const MatrixType a = MatrixType::Random(num_cols, num_cols);
@@ -97,24 +109,70 @@
   }
 }
 
-namespace {
-
-// NOTE: preprocessor directives in a macro are not standard conforming
-decltype(auto) MakeValues() {
-  return ::testing::Values(EIGEN
+INSTANTIATE_TEST_SUITE_P(
+    EigenCholesky,
+    DenseCholeskyTest,
+    ::testing::Combine(::testing::Values(EIGEN),
+                       ::testing::Values(kFullPrecision)),
+    ParamInfoToString);
 #ifndef CERES_NO_LAPACK
-                           ,
-                           LAPACK
+INSTANTIATE_TEST_SUITE_P(
+    LapackCholesky,
+    DenseCholeskyTest,
+    ::testing::Combine(::testing::Values(LAPACK),
+                       ::testing::Values(kFullPrecision)),
+    ParamInfoToString);
 #endif
 #ifndef CERES_NO_CUDA
-                           ,
-                           CUDA
+INSTANTIATE_TEST_SUITE_P(
+    CudaCholesky,
+    DenseCholeskyTest,
+    ::testing::Combine(::testing::Values(CUDA),
+                       ::testing::Values(kMixedPrecision,
+                                         kFullPrecision)),
+    ParamInfoToString);
 #endif
-  );
+
+TEST(DenseCholesky, ValidMixedPrecisionOptions) {
+#ifndef CERES_NO_CUDA
+  {
+    // Dense Cholesky with CUDA: okay, supported.
+    ContextImpl context;
+    LinearSolver::Options options;
+    options.dense_linear_algebra_library_type = CUDA;
+    options.use_mixed_precision_solves = true;
+    options.context = &context;
+    std::unique_ptr<DenseCholesky> dense_cholesky =
+        DenseCholesky::Create(options);
+    EXPECT_NE(dense_cholesky, nullptr);
+  }
+#endif
 }
 
-}  // namespace
+TEST(DenseCholesky, InvalidMixedPrecisionOptions) {
+  {
+    // Dense Cholesky with Eigen: not supported
+    ContextImpl context;
+    LinearSolver::Options options;
+    options.dense_linear_algebra_library_type = EIGEN;
+    options.use_mixed_precision_solves = true;
+    options.context = &context;
+    std::unique_ptr<DenseCholesky> dense_cholesky =
+        DenseCholesky::Create(options);
+    EXPECT_EQ(dense_cholesky, nullptr);
+  }
 
-INSTANTIATE_TEST_SUITE_P(_, DenseCholeskyTest, MakeValues(), ParamInfoToString);
+  {
+    // Dense Cholesky with Lapack: not supported
+    ContextImpl context;
+    LinearSolver::Options options;
+    options.dense_linear_algebra_library_type = LAPACK;
+    options.use_mixed_precision_solves = true;
+    options.context = &context;
+    std::unique_ptr<DenseCholesky> dense_cholesky =
+        DenseCholesky::Create(options);
+    EXPECT_EQ(dense_cholesky, nullptr);
+  }
+}
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index 9503843..93551f7 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -118,6 +118,48 @@
            internal::EigenSparse::IsNestedDissectionAvailable()));
 }
 
+bool MixedPrecisionOptionIsValid(const Solver::Options& options,
+                                 string* error) {
+  if (options.use_mixed_precision_solves) {
+    if ((options.linear_solver_type == DENSE_NORMAL_CHOLESKY ||
+        options.linear_solver_type == DENSE_SCHUR) &&
+        options.dense_linear_algebra_library_type == CUDA) {
+      // Mixed precision with CUDA and dense Cholesky variant: okay.
+      return true;
+    }
+    if ((options.linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        options.linear_solver_type == SPARSE_SCHUR) &&
+        (options.sparse_linear_algebra_library_type == EIGEN_SPARSE ||
+        options.sparse_linear_algebra_library_type == ACCELERATE_SPARSE)) {
+      // Mixed precision with any Eigen or Accelerate Cholesky variant: okay.
+      return true;
+    }
+    // No other mixed precision variants are supported.
+    if (options.linear_solver_type == DENSE_NORMAL_CHOLESKY ||
+        options.linear_solver_type == DENSE_SCHUR) {
+      *error = StringPrintf(
+          "use_mixed_precision_solves with %s is only supported with "
+          "CUDA as the dense_linear_algebra_library_type.",
+          LinearSolverTypeToString(options.linear_solver_type));
+      return false;
+    }
+    if ((options.linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        options.linear_solver_type == SPARSE_SCHUR) &&
+        options.sparse_linear_algebra_library_type == SUITE_SPARSE) {
+      *error =  StringPrintf(
+          "use_mixed_precision_solves with %s is not supported with "
+          "SUITE_SPARSE as the sparse_linear_algebra_library_type.",
+          LinearSolverTypeToString(options.linear_solver_type));
+      return false;
+    }
+    *error = StringPrintf(
+          "use_mixed_precision_solves with %s is not supported.",
+          LinearSolverTypeToString(options.linear_solver_type));
+    return false;
+  }
+  return true;
+}
+
 bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) {
   OPTION_GT(initial_trust_region_radius, 0.0);
   OPTION_GT(min_trust_region_radius, 0.0);
@@ -228,6 +270,10 @@
     }
   }
 
+  if (!MixedPrecisionOptionIsValid(options, error)) {
+    return false;
+  }
+
   if (options.trust_region_strategy_type == DOGLEG) {
     if (options.linear_solver_type == ITERATIVE_SCHUR ||
         options.linear_solver_type == CGNR) {
@@ -372,7 +418,7 @@
                                  &(summary->inner_iteration_ordering_given));
 
   // clang-format off
-  summary->dense_linear_algebra_library_type  = options.dense_linear_algebra_library_type;  
+  summary->dense_linear_algebra_library_type  = options.dense_linear_algebra_library_type;
   summary->dogleg_type                        = options.dogleg_type;
   summary->inner_iteration_time_in_seconds    = 0.0;
   summary->num_line_search_steps              = 0;
@@ -381,18 +427,18 @@
   summary->line_search_polynomial_minimization_time_in_seconds = 0.0;
   summary->line_search_total_time_in_seconds  = 0.0;
   summary->inner_iterations_given             = options.use_inner_iterations;
-  summary->line_search_direction_type         = options.line_search_direction_type;       
-  summary->line_search_interpolation_type     = options.line_search_interpolation_type;   
+  summary->line_search_direction_type         = options.line_search_direction_type;
+  summary->line_search_interpolation_type     = options.line_search_interpolation_type;
   summary->line_search_type                   = options.line_search_type;
   summary->linear_solver_type_given           = options.linear_solver_type;
   summary->max_lbfgs_rank                     = options.max_lbfgs_rank;
   summary->minimizer_type                     = options.minimizer_type;
-  summary->nonlinear_conjugate_gradient_type  = options.nonlinear_conjugate_gradient_type; 
+  summary->nonlinear_conjugate_gradient_type  = options.nonlinear_conjugate_gradient_type;
   summary->num_threads_given                  = options.num_threads;
   summary->preconditioner_type_given          = options.preconditioner_type;
   summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type;
-  summary->linear_solver_ordering_type        = options.linear_solver_ordering_type;       
-  summary->trust_region_strategy_type         = options.trust_region_strategy_type;        
+  summary->linear_solver_ordering_type        = options.linear_solver_ordering_type;
+  summary->trust_region_strategy_type         = options.trust_region_strategy_type;
   summary->visibility_clustering_type         = options.visibility_clustering_type;
   // clang-format on
 }
@@ -408,8 +454,9 @@
                                  &(summary->inner_iteration_ordering_used));
 
   // clang-format off
-  summary->inner_iterations_used          = pp.inner_iteration_minimizer != nullptr;   
+  summary->inner_iterations_used          = pp.inner_iteration_minimizer != nullptr;
   summary->linear_solver_type_used        = pp.linear_solver_options.type;
+  summary->mixed_precision_solves_used    = pp.options.use_mixed_precision_solves;
   summary->num_threads_used               = pp.options.num_threads;
   summary->preconditioner_type_used       = pp.options.preconditioner_type;
   // clang-format on
@@ -699,10 +746,13 @@
     if (linear_solver_type_used == DENSE_NORMAL_CHOLESKY ||
         linear_solver_type_used == DENSE_SCHUR ||
         linear_solver_type_used == DENSE_QR) {
+      const char* mixed_precision_suffix =
+          (mixed_precision_solves_used ? "(Mixed Precision)" : "");
       StringAppendF(&report,
-                    "\nDense linear algebra library  %15s\n",
+                    "\nDense linear algebra library  %15s %s\n",
                     DenseLinearAlgebraLibraryTypeToString(
-                        dense_linear_algebra_library_type));
+                        dense_linear_algebra_library_type),
+                    mixed_precision_suffix);
     }
 
     StringAppendF(&report,
@@ -723,12 +773,15 @@
         (linear_solver_type_used == ITERATIVE_SCHUR &&
          (preconditioner_type_used == CLUSTER_JACOBI ||
           preconditioner_type_used == CLUSTER_TRIDIAGONAL))) {
+      const char* mixed_precision_suffix =
+          (mixed_precision_solves_used ? "(Mixed Precision)" : "");
       StringAppendF(
           &report,
-          "\nSparse linear algebra library %15s + %s\n",
+          "\nSparse linear algebra library %15s + %s %s\n",
           SparseLinearAlgebraLibraryTypeToString(
               sparse_linear_algebra_library_type),
-          LinearSolverOrderingTypeToString(linear_solver_ordering_type));
+          LinearSolverOrderingTypeToString(linear_solver_ordering_type),
+          mixed_precision_suffix);
     }
 
     StringAppendF(&report, "\n");