diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml
index e4b124d..be4bad1 100644
--- a/.github/workflows/linux.yml
+++ b/.github/workflows/linux.yml
@@ -49,20 +49,23 @@
             libsuitesparse-dev \
             ninja-build \
             wget
-      - name: Setup CUDA toolkit (system repositories)
+      # nvidia cuda toolkit shipped with 20.04 LTS does not support stream-ordered allocations
+      - name: Setup CUDA toolkit repositories (20.04)
         if: matrix.gpu == 'cuda' && matrix.os == 'ubuntu:20.04'
         run: |
-          apt-get install -y \
-            nvidia-cuda-dev \
-            nvidia-cuda-toolkit
+          wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-keyring_1.0-1_all.deb
+          dpkg -i cuda-keyring_1.0-1_all.deb
       # nvidia cuda toolkit + gcc combo shipped with 22.04LTS is broken
       # and is not able to compile code that uses thrust
       # https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=1006962
-      - name: Setup CUDA toolkit (nvidia repositories)
+      - name: Setup CUDA toolkit repositories (22.04)
         if: matrix.gpu == 'cuda' && matrix.os == 'ubuntu:22.04'
         run: |
           wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2204/x86_64/cuda-keyring_1.0-1_all.deb
           dpkg -i cuda-keyring_1.0-1_all.deb
+      - name: Setup CUDA toolkit
+        if: matrix.gpu == 'cuda'
+        run: |
           apt-get update
           apt-get install -y cuda
           echo "CUDACXX=/usr/local/cuda/bin/nvcc" >> $GITHUB_ENV
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 9cb99b4..b67b10e 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -140,9 +140,12 @@
         COMMAND ${CUDAToolkit_BIN_DIR}/cuda-memcheck --leak-check full
             $<TARGET_FILE:cuda_dense_cholesky_test>)
   endif (BUILD_TESTING AND GFLAGS)
-  set_source_files_properties(cuda_kernels.cu.cc PROPERTIES LANGUAGE CUDA)
-  add_library(ceres_cuda_kernels STATIC cuda_kernels.cu.cc)
+  set_source_files_properties(cuda_kernels_vector_ops.cu.cc PROPERTIES LANGUAGE CUDA)
+  set_source_files_properties(cuda_kernels_bsm_to_crs.cu.cc PROPERTIES LANGUAGE CUDA)
+  add_library(ceres_cuda_kernels STATIC cuda_kernels_vector_ops.cu.cc cuda_kernels_bsm_to_crs.cu.cc)
   target_compile_features(ceres_cuda_kernels PRIVATE cxx_std_14)
+  # Enable __host__ / __device__ annotations in lambda declarations
+  target_compile_options(ceres_cuda_kernels PRIVATE --extended-lambda)
   target_include_directories(ceres_cuda_kernels PRIVATE  ${Ceres_SOURCE_DIR}/internal  ${Ceres_BINARY_DIR}/${CMAKE_INSTALL_INCLUDEDIR})
   list(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ceres_cuda_kernels)
 endif (USE_CUDA)
@@ -479,7 +482,7 @@
   ceres_test(cuda_block_structure)
   ceres_test(cuda_dense_cholesky)
   ceres_test(cuda_dense_qr)
-  ceres_test(cuda_kernels)
+  ceres_test(cuda_kernels_vector_ops)
   ceres_test(cuda_sparse_matrix)
   ceres_test(cuda_streamed_buffer)
   ceres_test(cuda_vector)
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index ab1d746..01dbfbe 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -709,7 +709,9 @@
 }
 
 std::unique_ptr<BlockSparseMatrix> BlockSparseMatrix::CreateRandomMatrix(
-    const BlockSparseMatrix::RandomMatrixOptions& options, std::mt19937& prng) {
+    const BlockSparseMatrix::RandomMatrixOptions& options,
+    std::mt19937& prng,
+    bool use_page_locked_memory) {
   CHECK_GT(options.num_row_blocks, 0);
   CHECK_GT(options.min_row_block_size, 0);
   CHECK_GT(options.max_row_block_size, 0);
@@ -766,7 +768,8 @@
     }
   }
 
-  auto matrix = std::make_unique<BlockSparseMatrix>(bs.release());
+  auto matrix =
+      std::make_unique<BlockSparseMatrix>(bs.release(), use_page_locked_memory);
   double* values = matrix->mutable_values();
   std::normal_distribution<double> standard_normal_distribution;
   std::generate_n(
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index 0d99e15..cd13b06 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -157,7 +157,9 @@
   // distributed and whose structure is determined by
   // RandomMatrixOptions.
   static std::unique_ptr<BlockSparseMatrix> CreateRandomMatrix(
-      const RandomMatrixOptions& options, std::mt19937& prng);
+      const RandomMatrixOptions& options,
+      std::mt19937& prng,
+      bool use_page_locked_memory = false);
 
  private:
   double* AllocateValues(int size);
diff --git a/internal/ceres/cuda_block_sparse_crs_view.cc b/internal/ceres/cuda_block_sparse_crs_view.cc
index a03d01c..bb267da 100644
--- a/internal/ceres/cuda_block_sparse_crs_view.cc
+++ b/internal/ceres/cuda_block_sparse_crs_view.cc
@@ -32,41 +32,69 @@
 
 #ifndef CERES_NO_CUDA
 
-#include "ceres/cuda_block_structure.h"
-#include "ceres/cuda_kernels.h"
+#include "ceres/cuda_kernels_bsm_to_crs.h"
 
 namespace ceres::internal {
 
 CudaBlockSparseCRSView::CudaBlockSparseCRSView(const BlockSparseMatrix& bsm,
                                                ContextImpl* context)
-    : permutation_(context, bsm.num_nonzeros()),
-      streamed_buffer_(context, kMaxTemporaryArraySize) {
-  CudaBlockSparseStructure block_structure(*bsm.block_structure(), context);
+    : context_(context) {
+  block_structure_ = std::make_unique<CudaBlockSparseStructure>(
+      *bsm.block_structure(), context);
+  // Only block-sparse matrices with sequential layout of cells are supported
+  CHECK(block_structure_->sequential_layout());
   crs_matrix_ = std::make_unique<CudaSparseMatrix>(
       bsm.num_rows(), bsm.num_cols(), bsm.num_nonzeros(), context);
-  CudaBuffer<int> temp_rows(context, bsm.num_rows());
-  FillCRSStructure(block_structure.num_row_blocks(),
+  FillCRSStructure(block_structure_->num_row_blocks(),
                    bsm.num_rows(),
-                   block_structure.row_block_offsets(),
-                   block_structure.cells(),
-                   block_structure.row_blocks(),
-                   block_structure.col_blocks(),
+                   block_structure_->first_cell_in_row_block(),
+                   block_structure_->cells(),
+                   block_structure_->row_blocks(),
+                   block_structure_->col_blocks(),
                    crs_matrix_->mutable_rows(),
                    crs_matrix_->mutable_cols(),
-                   temp_rows.data(),
-                   permutation_.data(),
                    context->DefaultStream());
+  is_crs_compatible_ = block_structure_->IsCrsCompatible();
+  // if matrix is crs-compatible - we can drop block-structure and don't need
+  // streamed_buffer_
+  if (is_crs_compatible_) {
+    VLOG(3) << "Block-sparse matrix is compatible with CRS, discarding "
+               "block-structure";
+    block_structure_ = nullptr;
+  } else {
+    streamed_buffer_ = std::make_unique<CudaStreamedBuffer<double>>(
+        context_, kMaxTemporaryArraySize);
+  }
   UpdateValues(bsm);
 }
+
 void CudaBlockSparseCRSView::UpdateValues(const BlockSparseMatrix& bsm) {
-  streamed_buffer_.CopyToGpu(
+  if (is_crs_compatible_) {
+    // Values of CRS-compatible matrices can be copied as-is
+    CHECK_EQ(cudaSuccess,
+             cudaMemcpyAsync(crs_matrix_->mutable_values(),
+                             bsm.values(),
+                             bsm.num_nonzeros() * sizeof(double),
+                             cudaMemcpyHostToDevice,
+                             context_->DefaultStream()));
+    return;
+  }
+  streamed_buffer_->CopyToGpu(
       bsm.values(),
       bsm.num_nonzeros(),
-      [permutation = permutation_.data(),
-       values_to = crs_matrix_->mutable_values()](
+      [bs = block_structure_.get(), crs = crs_matrix_.get()](
           const double* values, int num_values, int offset, auto stream) {
-        PermuteValues(
-            offset, num_values, permutation, values, values_to, stream);
+        PermuteToCRS(offset,
+                     num_values,
+                     bs->num_row_blocks(),
+                     bs->first_cell_in_row_block(),
+                     bs->cells(),
+                     bs->row_blocks(),
+                     bs->col_blocks(),
+                     crs->rows(),
+                     values,
+                     crs->mutable_values(),
+                     stream);
       });
 }
 
diff --git a/internal/ceres/cuda_block_sparse_crs_view.h b/internal/ceres/cuda_block_sparse_crs_view.h
index 3ea8498..2ae8721 100644
--- a/internal/ceres/cuda_block_sparse_crs_view.h
+++ b/internal/ceres/cuda_block_sparse_crs_view.h
@@ -39,6 +39,7 @@
 #include <memory>
 
 #include "ceres/block_sparse_matrix.h"
+#include "ceres/cuda_block_structure.h"
 #include "ceres/cuda_buffer.h"
 #include "ceres/cuda_sparse_matrix.h"
 #include "ceres/cuda_streamed_buffer.h"
@@ -49,42 +50,23 @@
 // following operations in order to compute products of block-sparse matrices
 // and dense vectors on gpu:
 //  - Once per block-sparse structure update:
-//    - Compute CRS structure from block-sparse structure
-//    - Compute permutation from block-sparse values to CRS values
+//    - Compute CRS structure from block-sparse structure and check if values of
+//    block-sparse matrix would have the same order as values of CRS matrix
 //  - Once per block-sparse values update:
 //    - Update values in CRS matrix with values of block-sparse matrix
 //
-// Since there are no constraints on positions of cells in value array of
-// block-sparse matrix, a permutation from block-sparse values to CRS
-// values is stored explicitly.
+// Only block-sparse matrices with sequential order of cells are supported.
 //
-// Example: given matrix with the following block-structure
-//  [ 1 2 | | 6 7 ]
-//  [ 3 4 | | 8 9 ]
-//  [-----+-+-----]
-//  [     |5|     ]
-// with values stored as values_block_sparse = [1, 2, 3, 4, 5, 6, 7, 8, 9],
-// permutation from block-sparse to CRS is p = [0, 1, 4, 5, 8, 2, 3, 6, 7];
-// and i-th block-sparse value has index p[i] in CRS values array
-//
-// This allows to avoid storing both CRS and block-sparse values in GPU memory.
-// Instead, block-sparse values are transferred to gpu memory as a disjoint set
-// of small continuous segments with simultaneous permutation of the values into
-// correct order
+// UpdateValues method updates values:
+//  - In a single host-to-device copy for matrices with CRS-compatible value
+//  layout
+//  - Simultaneously transferring and permuting values using CudaStreamedBuffer
+//  otherwise
 class CERES_NO_EXPORT CudaBlockSparseCRSView {
  public:
-  // Initializes internal CRS matrix and permutation from block-sparse to CRS
-  // values. The following objects are stored in gpu memory for the whole
-  // lifetime of the object
-  //  - crs_matrix_: CRS matrix
-  //  - permutation_: permutation from block-sparse to CRS value order
-  //  (num_nonzeros integer values)
-  //  - streamed_buffer_: helper for value updating
-  // The following objects are created temporarily during construction:
-  //  - CudaBlockSparseStructure: block-sparse structure of block-sparse matrix
-  //  - num_rows integer values: row to row-block map
-  //  If copy_values flag is set to false, only structure of block-sparse matrix
-  //  bsm is captured, and values are left uninitialized
+  // Initializes internal CRS matrix using structure and values of block-sparse
+  // matrix For block-sparse matrices that have value layout different from CRS
+  // block-sparse structure will be stored/
   CudaBlockSparseCRSView(const BlockSparseMatrix& bsm, ContextImpl* context);
 
   const CudaSparseMatrix* crs_matrix() const { return crs_matrix_.get(); }
@@ -95,16 +77,21 @@
   // used for construction.
   void UpdateValues(const BlockSparseMatrix& bsm);
 
+  // Returns true if block-sparse matrix had CRS-compatible value layout
+  bool IsCrsCompatible() const { return is_crs_compatible_; }
+
  private:
   // Value permutation kernel performs a single element-wise operation per
   // thread, thus performing permutation in blocks of 8 megabytes of
   // block-sparse  values seems reasonable
   static constexpr int kMaxTemporaryArraySize = 1 * 1024 * 1024;
   std::unique_ptr<CudaSparseMatrix> crs_matrix_;
-  // Permutation from block-sparse to CRS value order.
-  // permutation_[i] = index of i-th block-sparse value in CRS values
-  CudaBuffer<int> permutation_;
-  CudaStreamedBuffer<double> streamed_buffer_;
+  // Only created if block-sparse matrix has non-CRS value layout
+  std::unique_ptr<CudaStreamedBuffer<double>> streamed_buffer_;
+  // Only stored if block-sparse matrix has non-CRS value layout
+  std::unique_ptr<CudaBlockSparseStructure> block_structure_;
+  bool is_crs_compatible_;
+  ContextImpl* context_;
 };
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/cuda_block_sparse_crs_view_test.cc b/internal/ceres/cuda_block_sparse_crs_view_test.cc
index f8d1cdf..e80d026 100644
--- a/internal/ceres/cuda_block_sparse_crs_view_test.cc
+++ b/internal/ceres/cuda_block_sparse_crs_view_test.cc
@@ -52,48 +52,110 @@
     options.max_col_block_size = 10;
     options.block_density = 0.2;
     std::mt19937 rng;
-    A_ = BlockSparseMatrix::CreateRandomMatrix(options, rng);
-    std::iota(
-        A_->mutable_values(), A_->mutable_values() + A_->num_nonzeros(), 1);
+
+    // Block-sparse matrix with order of values different from CRS
+    block_sparse_non_crs_compatible_ =
+        BlockSparseMatrix::CreateRandomMatrix(options, rng, true);
+    std::iota(block_sparse_non_crs_compatible_->mutable_values(),
+              block_sparse_non_crs_compatible_->mutable_values() +
+                  block_sparse_non_crs_compatible_->num_nonzeros(),
+              1);
+
+    options.max_row_block_size = 1;
+    // Block-sparse matrix with CRS order of values (row-blocks are rows)
+    block_sparse_crs_compatible_rows_ =
+        BlockSparseMatrix::CreateRandomMatrix(options, rng, true);
+    std::iota(block_sparse_crs_compatible_rows_->mutable_values(),
+              block_sparse_crs_compatible_rows_->mutable_values() +
+                  block_sparse_crs_compatible_rows_->num_nonzeros(),
+              1);
+    // Block-sparse matrix with CRS order of values (single cell per row-block)
+    auto bs = std::make_unique<CompressedRowBlockStructure>(
+        *block_sparse_non_crs_compatible_->block_structure());
+
+    int num_nonzeros = 0;
+    for (auto& r : bs->rows) {
+      const int num_cells = r.cells.size();
+      if (num_cells > 1) {
+        std::uniform_int_distribution<int> uniform_cell(0, num_cells - 1);
+        const int selected_cell = uniform_cell(rng);
+        std::swap(r.cells[0], r.cells[selected_cell]);
+        r.cells.resize(1);
+      }
+      const int row_block_size = r.block.size;
+      for (auto& c : r.cells) {
+        c.position = num_nonzeros;
+        const int col_block_size = bs->cols[c.block_id].size;
+        num_nonzeros += col_block_size * row_block_size;
+      }
+    }
+    block_sparse_crs_compatible_single_cell_ =
+        std::make_unique<BlockSparseMatrix>(bs.release());
+    std::iota(block_sparse_crs_compatible_single_cell_->mutable_values(),
+              block_sparse_crs_compatible_single_cell_->mutable_values() +
+                  block_sparse_crs_compatible_single_cell_->num_nonzeros(),
+              1);
   }
 
-  std::unique_ptr<BlockSparseMatrix> A_;
+  void Compare(const BlockSparseMatrix& bsm, CudaSparseMatrix& csm) {
+    ASSERT_EQ(csm.num_cols(), bsm.num_cols());
+    ASSERT_EQ(csm.num_rows(), bsm.num_rows());
+    ASSERT_EQ(csm.num_nonzeros(), bsm.num_nonzeros());
+    const int num_rows = bsm.num_rows();
+    const int num_cols = bsm.num_cols();
+    Vector x(num_cols);
+    Vector y(num_rows);
+    CudaVector x_cuda(&context_, num_cols);
+    CudaVector y_cuda(&context_, num_rows);
+    Vector y_cuda_host(num_rows);
+
+    for (int i = 0; i < num_cols; ++i) {
+      x.setZero();
+      y.setZero();
+      y_cuda.SetZero();
+      x[i] = 1.;
+      x_cuda.CopyFromCpu(x);
+      csm.RightMultiplyAndAccumulate(x_cuda, &y_cuda);
+      bsm.RightMultiplyAndAccumulate(
+          x.data(), y.data(), &context_, std::thread::hardware_concurrency());
+      y_cuda.CopyTo(&y_cuda_host);
+      // There will be up to 1 non-zero product per row, thus we expect an exact
+      // match on 32-bit integer indices
+      EXPECT_EQ((y - y_cuda_host).squaredNorm(), 0.);
+    }
+  }
+
+  std::unique_ptr<BlockSparseMatrix> block_sparse_non_crs_compatible_;
+  std::unique_ptr<BlockSparseMatrix> block_sparse_crs_compatible_rows_;
+  std::unique_ptr<BlockSparseMatrix> block_sparse_crs_compatible_single_cell_;
   ContextImpl context_;
 };
 
-TEST_F(CudaBlockSparseCRSViewTest, CreateUpdateValues) {
-  auto view = CudaBlockSparseCRSView(*A_, &context_);
+TEST_F(CudaBlockSparseCRSViewTest, CreateUpdateValuesNonCompatible) {
+  auto view =
+      CudaBlockSparseCRSView(*block_sparse_non_crs_compatible_, &context_);
+  ASSERT_EQ(view.IsCrsCompatible(), false);
 
-  // TODO: we definitely would like to use crs_matrix() here, but
-  // CudaSparseMatrix::RightMultiplyAndAccumulate is defined non-const because
-  // it might allocate additional storage by request of cuSPARSE
   auto matrix = view.mutable_crs_matrix();
-  ASSERT_EQ(matrix->num_cols(), A_->num_cols());
-  ASSERT_EQ(matrix->num_rows(), A_->num_rows());
-  ASSERT_EQ(matrix->num_nonzeros(), A_->num_nonzeros());
+  Compare(*block_sparse_non_crs_compatible_, *matrix);
+}
 
-  const int num_rows = A_->num_rows();
-  const int num_cols = A_->num_cols();
-  Vector x(num_cols);
-  Vector y(num_rows);
-  CudaVector x_cuda(&context_, num_cols);
-  CudaVector y_cuda(&context_, num_rows);
-  Vector y_cuda_host(num_rows);
+TEST_F(CudaBlockSparseCRSViewTest, CreateUpdateValuesCompatibleRows) {
+  auto view =
+      CudaBlockSparseCRSView(*block_sparse_crs_compatible_rows_, &context_);
+  ASSERT_EQ(view.IsCrsCompatible(), true);
 
-  for (int i = 0; i < num_cols; ++i) {
-    x.setZero();
-    y.setZero();
-    y_cuda.SetZero();
-    x[i] = 1.;
-    x_cuda.CopyFromCpu(x);
-    A_->RightMultiplyAndAccumulate(
-        x.data(), y.data(), &context_, std::thread::hardware_concurrency());
-    matrix->RightMultiplyAndAccumulate(x_cuda, &y_cuda);
-    y_cuda.CopyTo(&y_cuda_host);
-    // There will be up to 1 non-zero product per row, thus we expect an exact
-    // match on 32-bit integer indices
-    EXPECT_EQ((y - y_cuda_host).squaredNorm(), 0.);
-  }
+  auto matrix = view.mutable_crs_matrix();
+  Compare(*block_sparse_crs_compatible_rows_, *matrix);
+}
+
+TEST_F(CudaBlockSparseCRSViewTest, CreateUpdateValuesCompatibleSingleCell) {
+  auto view = CudaBlockSparseCRSView(*block_sparse_crs_compatible_single_cell_,
+                                     &context_);
+  ASSERT_EQ(view.IsCrsCompatible(), true);
+
+  auto matrix = view.mutable_crs_matrix();
+  Compare(*block_sparse_crs_compatible_single_cell_, *matrix);
 }
 }  // namespace ceres::internal
 
diff --git a/internal/ceres/cuda_block_structure.cc b/internal/ceres/cuda_block_structure.cc
index 8e549e1..5817cd7 100644
--- a/internal/ceres/cuda_block_structure.cc
+++ b/internal/ceres/cuda_block_structure.cc
@@ -46,7 +46,7 @@
 
 CudaBlockSparseStructure::CudaBlockSparseStructure(
     const CompressedRowBlockStructure& block_structure, ContextImpl* context)
-    : row_block_offsets_(context),
+    : first_cell_in_row_block_(context),
       cells_(context),
       row_blocks_(context),
       col_blocks_(context) {
@@ -56,44 +56,56 @@
   const auto& col_blocks = block_structure.cols;
 
   // Row block offset is an index of the first cell corresponding to row block
-  std::vector<int> row_block_offsets;
+  std::vector<int> first_cell_in_row_block;
   // Flat array of all cells from all row-blocks
   std::vector<Cell> cells;
 
+  int f_values_offset = 0;
+  is_crs_compatible_ = true;
   num_row_blocks_ = block_structure.rows.size();
   num_col_blocks_ = col_blocks.size();
 
   row_blocks.reserve(num_row_blocks_);
-  row_block_offsets.reserve(num_row_blocks_ + 1);
+  first_cell_in_row_block.reserve(num_row_blocks_ + 1);
   num_nonzeros_ = 0;
-  num_cells_ = 0;
+  sequential_layout_ = true;
   for (const auto& r : block_structure.rows) {
     const int row_block_size = r.block.size;
+    if (r.cells.size() > 1 && row_block_size > 1) {
+      is_crs_compatible_ = false;
+    }
     row_blocks.emplace_back(r.block);
-    row_block_offsets.push_back(num_cells_);
+    first_cell_in_row_block.push_back(cells.size());
     for (const auto& c : r.cells) {
-      cells.emplace_back(c);
       const int col_block_size = col_blocks[c.block_id].size;
-      num_nonzeros_ += col_block_size * row_block_size;
-      ++num_cells_;
+      const int cell_size = col_block_size * row_block_size;
+      cells.push_back(c);
+      sequential_layout_ &= c.position == num_nonzeros_;
+      num_nonzeros_ += cell_size;
     }
   }
-  row_block_offsets.push_back(num_cells_);
+  first_cell_in_row_block.push_back(cells.size());
+  num_cells_ = cells.size();
 
   num_rows_ = Dimension(row_blocks);
   num_cols_ = Dimension(col_blocks);
 
+  is_crs_compatible_ &= sequential_layout_;
+
   if (VLOG_IS_ON(3)) {
-    const size_t row_block_offsets_size =
-        row_block_offsets.size() * sizeof(int);
+    const size_t first_cell_in_row_block_size =
+        first_cell_in_row_block.size() * sizeof(int);
     const size_t cells_size = cells.size() * sizeof(Cell);
     const size_t row_blocks_size = row_blocks.size() * sizeof(Block);
     const size_t col_blocks_size = col_blocks.size() * sizeof(Block);
-    const size_t total_size =
-        row_block_offsets_size + cells_size + col_blocks_size + row_blocks_size;
+    const size_t total_size = first_cell_in_row_block_size + cells_size +
+                              col_blocks_size + row_blocks_size;
+    const double ratio =
+        (100. * total_size) / (num_nonzeros_ * (sizeof(int) + sizeof(double)) +
+                               num_rows_ * sizeof(int));
     VLOG(3) << "\nCudaBlockSparseStructure:\n"
                "\tRow block offsets: "
-            << row_block_offsets_size
+            << first_cell_in_row_block_size
             << " bytes\n"
                "\tColumn blocks: "
             << col_blocks_size
@@ -102,13 +114,11 @@
             << row_blocks_size
             << " bytes\n"
                "\tCells: "
-            << cells_size
-            << " bytes\n"
-               "\tTotal: "
-            << total_size << " bytes of GPU memory";
+            << cells_size << " bytes\n\tTotal: " << total_size
+            << " bytes of GPU memory (" << ratio << "% of CRS matrix size)";
   }
 
-  row_block_offsets_.CopyFromCpuVector(row_block_offsets);
+  first_cell_in_row_block_.CopyFromCpuVector(first_cell_in_row_block);
   cells_.CopyFromCpuVector(cells);
   row_blocks_.CopyFromCpuVector(row_blocks);
   col_blocks_.CopyFromCpuVector(col_blocks);
diff --git a/internal/ceres/cuda_block_structure.h b/internal/ceres/cuda_block_structure.h
index a41d0b3..28e36b2 100644
--- a/internal/ceres/cuda_block_structure.h
+++ b/internal/ceres/cuda_block_structure.h
@@ -44,12 +44,11 @@
 // This class stores a read-only block-sparse structure in gpu memory.
 // Invariants are the same as those of CompressedRowBlockStructure.
 // In order to simplify allocation and copying data to gpu, cells from all
-// row-blocks are stored in a single array sequentially. Additional array
-// row_block_offsets of size num_row_blocks + 1 allows to identify range of
-// cells corresponding to a row-block.
-// Cells corresponding to i-th row-block are stored in sub-array
-// cells[row_block_offsets[i]; ... row_block_offsets[i + 1] - 1], and their
-// order is preserved.
+// row-blocks are stored in a single array sequentially. Array
+// first_cell_in_row_block of size num_row_blocks + 1 allows to identify range
+// of cells corresponding to a row-block. Cells corresponding to i-th row-block
+// are stored in sub-array cells[first_cell_in_row_block[i]; ...
+// first_cell_in_row_block[i + 1] - 1], and their order is preserved.
 class CERES_NO_EXPORT CudaBlockSparseStructure {
  public:
   // CompressedRowBlockStructure is contains a vector of CompressedLists, with
@@ -65,9 +64,20 @@
   int num_row_blocks() const { return num_row_blocks_; }
   int num_col_blocks() const { return num_col_blocks_; }
 
+  // Returns true if values from block-sparse matrix can be copied to CRS matrix
+  // as-is. This is possible if each row-block is stored in CRS order:
+  //  - Row-block consists of a single row
+  //  - Row-block contains a single cell
+  bool IsCrsCompatible() const { return is_crs_compatible_; }
+  // Returns true if block-sparse structure corresponds to block-sparse matrix
+  // with sequential cell positions
+  bool sequential_layout() const { return sequential_layout_; }
+
   // Device pointer to array of num_row_blocks + 1 indices of the first cell of
   // row block
-  const int* row_block_offsets() const { return row_block_offsets_.data(); }
+  const int* first_cell_in_row_block() const {
+    return first_cell_in_row_block_.data();
+  }
   // Device pointer to array of num_cells cells, sorted by row-block
   const Cell* cells() const { return cells_.data(); }
   // Device pointer to array of row blocks
@@ -82,7 +92,9 @@
   int num_nonzeros_;
   int num_row_blocks_;
   int num_col_blocks_;
-  CudaBuffer<int> row_block_offsets_;
+  bool is_crs_compatible_;
+  bool sequential_layout_;
+  CudaBuffer<int> first_cell_in_row_block_;
   CudaBuffer<Cell> cells_;
   CudaBuffer<Block> row_blocks_;
   CudaBuffer<Block> col_blocks_;
diff --git a/internal/ceres/cuda_block_structure_test.cc b/internal/ceres/cuda_block_structure_test.cc
index 95cb26b..000b6bb 100644
--- a/internal/ceres/cuda_block_structure_test.cc
+++ b/internal/ceres/cuda_block_structure_test.cc
@@ -83,10 +83,10 @@
   }
   std::vector<int> GetRowBlockOffsets(
       const CudaBlockSparseStructure& structure) {
-    const auto& cuda_buffer = structure.row_block_offsets_;
-    std::vector<int> row_block_offsets(cuda_buffer.size());
-    cuda_buffer.CopyToCpu(row_block_offsets.data(), row_block_offsets.size());
-    return row_block_offsets;
+    const auto& cuda_buffer = structure.first_cell_in_row_block_;
+    std::vector<int> first_cell_in_row_block(cuda_buffer.size());
+    cuda_buffer.CopyToCpu(first_cell_in_row_block.data(), first_cell_in_row_block.size());
+    return first_cell_in_row_block;
   }
 
   std::unique_ptr<BlockSparseMatrix> A_;
@@ -114,19 +114,19 @@
   }
 
   std::vector<Cell> cells = GetCells(cuda_block_structure);
-  std::vector<int> row_block_offsets = GetRowBlockOffsets(cuda_block_structure);
+  std::vector<int> first_cell_in_row_block = GetRowBlockOffsets(cuda_block_structure);
   blocks = GetRowBlocks(cuda_block_structure);
 
   ASSERT_EQ(blocks.size(), num_row_blocks);
-  ASSERT_EQ(row_block_offsets.size(), num_row_blocks + 1);
-  ASSERT_EQ(row_block_offsets.back(), cells.size());
+  ASSERT_EQ(first_cell_in_row_block.size(), num_row_blocks + 1);
+  ASSERT_EQ(first_cell_in_row_block.back(), cells.size());
 
   for (int i = 0; i < num_row_blocks; ++i) {
     const int num_cells = block_structure->rows[i].cells.size();
     EXPECT_EQ(blocks[i].position, block_structure->rows[i].block.position);
     EXPECT_EQ(blocks[i].size, block_structure->rows[i].block.size);
-    const int first_cell = row_block_offsets[i];
-    const int last_cell = row_block_offsets[i + 1];
+    const int first_cell = first_cell_in_row_block[i];
+    const int last_cell = first_cell_in_row_block[i + 1];
     ASSERT_EQ(last_cell - first_cell, num_cells);
     for (int j = 0; j < num_cells; ++j) {
       EXPECT_EQ(cells[first_cell + j].block_id,
diff --git a/internal/ceres/cuda_kernels.cu.cc b/internal/ceres/cuda_kernels.cu.cc
deleted file mode 100644
index 3d93ae8..0000000
--- a/internal/ceres/cuda_kernels.cu.cc
+++ /dev/null
@@ -1,304 +0,0 @@
-// 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/cuda_kernels.h"
-
-#include <thrust/execution_policy.h>
-#include <thrust/scan.h>
-
-#include "block_structure.h"
-#include "cuda_runtime.h"
-
-namespace ceres {
-namespace 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;
-inline int NumBlocks(int size) {
-  return (size + kCudaBlockSize - 1) / kCudaBlockSize;
-}
-
-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 = NumBlocks(size);
-  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 = NumBlocks(size);
-  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 = NumBlocks(size);
-  SetZeroKernel<float><<<num_blocks, kCudaBlockSize, 0, stream>>>(output, size);
-}
-
-void CudaSetZeroFP64(double* output, const int size, cudaStream_t stream) {
-  const int num_blocks = NumBlocks(size);
-  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 = NumBlocks(size);
-  XPlusEqualsYKernel<float, double>
-      <<<num_blocks, kCudaBlockSize, 0, stream>>>(x, y, size);
-}
-
-__global__ void CudaDtDxpyKernel(double* __restrict__ y,
-                                 const double* D,
-                                 const double* __restrict__ x,
-                                 const int size) {
-  const int i = blockIdx.x * blockDim.x + threadIdx.x;
-  if (i < size) {
-    y[i] = y[i] + D[i] * D[i] * x[i];
-  }
-}
-
-void CudaDtDxpy(double* y,
-                const double* D,
-                const double* x,
-                const int size,
-                cudaStream_t stream) {
-  const int num_blocks = NumBlocks(size);
-  CudaDtDxpyKernel<<<num_blocks, kCudaBlockSize, 0, stream>>>(y, D, x, size);
-}
-
-// Fill row block id and nnz for each row using block-sparse structure
-// represented by a set of flat arrays.
-// Inputs:
-// - num_row_blocks: number of row-blocks in block-sparse structure
-// - row_block_offsets: index of the first cell of the row-block; size:
-// num_row_blocks + 1
-// - cells: cells of block-sparse structure as a continuous array
-// - row_blocks: row blocks of block-sparse structure stored sequentially
-// - col_blocks: column blocks of block-sparse structure stored sequentially
-// Outputs:
-// - rows: rows[i + 1] will contain number of non-zeros in i-th row, rows[0]
-// will be set to 0; rows are filled with a shift by one element in order
-// to obtain row-index array of CRS matrix with a inclusive scan afterwards
-// - row_block_ids: row_block_ids[i] will be set to index of row-block that
-// contains i-th row.
-// Computation is perform row-block-wise
-__global__ void RowBlockIdAndNNZ(int num_row_blocks,
-                                 const int* __restrict__ row_block_offsets,
-                                 const Cell* __restrict__ cells,
-                                 const Block* __restrict__ row_blocks,
-                                 const Block* __restrict__ col_blocks,
-                                 int* __restrict__ rows,
-                                 int* __restrict__ row_block_ids) {
-  const int row_block_id = blockIdx.x * blockDim.x + threadIdx.x;
-  if (row_block_id > num_row_blocks) {
-    // No synchronization is performed in this kernel, thus it is safe to return
-    return;
-  }
-  if (row_block_id == num_row_blocks) {
-    // one extra thread sets the first element
-    rows[0] = 0;
-    return;
-  }
-  const auto& row_block = row_blocks[row_block_id];
-  int row_nnz = 0;
-  const auto first_cell = cells + row_block_offsets[row_block_id];
-  const auto last_cell = cells + row_block_offsets[row_block_id + 1];
-  for (auto cell = first_cell; cell < last_cell; ++cell) {
-    row_nnz += col_blocks[cell->block_id].size;
-  }
-  const int first_row = row_block.position;
-  const int last_row = first_row + row_block.size;
-  for (int i = first_row; i < last_row; ++i) {
-    rows[i + 1] = row_nnz;
-    row_block_ids[i] = row_block_id;
-  }
-}
-
-// Row-wise creation of CRS structure
-// Inputs:
-// - num_rows: number of rows in matrix
-// - row_block_offsets: index of the first cell of the row-block; size:
-// num_row_blocks + 1
-// - cells: cells of block-sparse structure as a continuous array
-// - row_blocks: row blocks of block-sparse structure stored sequentially
-// - col_blocks: column blocks of block-sparse structure stored sequentially
-// - row_block_ids: index of row-block that corresponds to row
-// - rows: row-index array of CRS structure
-// Outputs:
-// - cols: column-index array of CRS structure
-// - permutation: permutation from block-sparse to crs order
-// Computaion is perform row-wise
-__global__ void ComputeColumnsAndPermutation(
-    int num_rows,
-    const int* __restrict__ row_block_offsets,
-    const Cell* __restrict__ cells,
-    const Block* __restrict__ row_blocks,
-    const Block* __restrict__ col_blocks,
-    const int* __restrict__ row_block_ids,
-    const int* __restrict__ rows,
-    int* __restrict__ cols,
-    int* __restrict__ permutation) {
-  const int row = blockIdx.x * blockDim.x + threadIdx.x;
-  if (row >= num_rows) {
-    // No synchronization is performed in this kernel, thus it is safe to return
-    return;
-  }
-  const int row_block_id = row_block_ids[row];
-  // zero-based index of row in row-block
-  const int row_in_block = row - row_blocks[row_block_id].position;
-  // position in crs matrix
-  int crs_position = rows[row];
-  const auto first_cell = cells + row_block_offsets[row_block_id];
-  const auto last_cell = cells + row_block_offsets[row_block_id + 1];
-  // For reach cell of row-block only current row is being filled
-  for (auto cell = first_cell; cell < last_cell; ++cell) {
-    const auto& col_block = col_blocks[cell->block_id];
-    const int col_block_size = col_block.size;
-    int column_idx = col_block.position;
-    int bs_position = cell->position + row_in_block * col_block_size;
-    // Fill permutation and column indices for each element of row_in_block
-    // row of current cell
-    for (int i = 0; i < col_block_size; ++i, ++crs_position) {
-      permutation[bs_position++] = crs_position;
-      cols[crs_position] = column_idx++;
-    }
-  }
-}
-
-void FillCRSStructure(const int num_row_blocks,
-                      const int num_rows,
-                      const int* row_block_offsets,
-                      const Cell* cells,
-                      const Block* row_blocks,
-                      const Block* col_blocks,
-                      int* rows,
-                      int* cols,
-                      int* row_block_ids,
-                      int* permutation,
-                      cudaStream_t stream) {
-  // Set number of non-zeros per row in rows array and row to row-block map in
-  // row_block_ids array
-  const int num_blocks_blockwise = NumBlocks(num_row_blocks + 1);
-  RowBlockIdAndNNZ<<<num_blocks_blockwise, kCudaBlockSize, 0, stream>>>(
-      num_row_blocks,
-      row_block_offsets,
-      cells,
-      row_blocks,
-      col_blocks,
-      rows,
-      row_block_ids);
-  // Finalize row-index array of CRS strucure by computing prefix sum
-  thrust::inclusive_scan(
-      thrust::cuda::par.on(stream), rows, rows + num_rows + 1, rows);
-
-  // Fill cols array of CRS structure and permutation from block-sparse to CRS
-  const int num_blocks_rowwise = NumBlocks(num_rows);
-  ComputeColumnsAndPermutation<<<num_blocks_rowwise,
-                                 kCudaBlockSize,
-                                 0,
-                                 stream>>>(num_rows,
-                                           row_block_offsets,
-                                           cells,
-                                           row_blocks,
-                                           col_blocks,
-                                           row_block_ids,
-                                           rows,
-                                           cols,
-                                           permutation);
-}
-
-// Updates values of CRS matrix with a continuous block of values of
-// block-sparse matrix. With permutation[i] being an index of i-th value of
-// block-sparse matrix in values of CRS matrix updating values is quite
-// efficient when performed element-wise (reads of permutation and values arrays
-// are coalesced when offset is divisable by 16)
-__global__ void PermuteValuesKernel(const int offset,
-                                    const int num_values,
-                                    const int* permutation,
-                                    const double* block_sparse_values,
-                                    double* crs_values) {
-  const int value_id = blockIdx.x * blockDim.x + threadIdx.x;
-  if (value_id < num_values) {
-    crs_values[permutation[offset + value_id]] = block_sparse_values[value_id];
-  }
-}
-
-void PermuteValues(const int offset,
-                   const int num_values,
-                   const int* permutation,
-                   const double* block_sparse_values,
-                   double* crs_values,
-                   cudaStream_t stream) {
-  const int num_blocks = NumBlocks(num_values);
-  PermuteValuesKernel<<<num_blocks, kCudaBlockSize, 0, stream>>>(
-      offset, num_values, permutation, block_sparse_values, crs_values);
-}
-
-}  // namespace internal
-}  // namespace ceres
diff --git a/internal/ceres/cuda_kernels_bsm_to_crs.cu.cc b/internal/ceres/cuda_kernels_bsm_to_crs.cu.cc
new file mode 100644
index 0000000..8be3553
--- /dev/null
+++ b/internal/ceres/cuda_kernels_bsm_to_crs.cu.cc
@@ -0,0 +1,279 @@
+// 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.
+//
+// Authors: dmitriy.korchemkin@gmail.com (Dmitriy Korchemkin)
+
+#include "ceres/cuda_kernels_bsm_to_crs.h"
+
+#include <cuda_runtime.h>
+#include <thrust/execution_policy.h>
+#include <thrust/scan.h>
+
+#include "ceres/block_structure.h"
+#include "ceres/cuda_kernels_utils.h"
+
+namespace ceres {
+namespace internal {
+
+// Fill row block id and nnz for each row using block-sparse structure
+// represented by a set of flat arrays.
+// Inputs:
+// - num_row_blocks: number of row-blocks in block-sparse structure
+// - first_cell_in_row_block: index of the first cell of the row-block; size:
+// num_row_blocks + 1
+// - cells: cells of block-sparse structure as a continuous array
+// - row_blocks: row blocks of block-sparse structure stored sequentially
+// - col_blocks: column blocks of block-sparse structure stored sequentially
+// Outputs:
+// - rows: rows[i + 1] will contain number of non-zeros in i-th row, rows[0]
+// will be set to 0; rows are filled with a shift by one element in order
+// to obtain row-index array of CRS matrix with a inclusive scan afterwards
+// - row_block_ids: row_block_ids[i] will be set to index of row-block that
+// contains i-th row.
+// Computation is perform row-block-wise
+__global__ void RowBlockIdAndNNZ(
+    int num_row_blocks,
+    const int* __restrict__ first_cell_in_row_block,
+    const Cell* __restrict__ cells,
+    const Block* __restrict__ row_blocks,
+    const Block* __restrict__ col_blocks,
+    int* __restrict__ rows,
+    int* __restrict__ row_block_ids) {
+  const int row_block_id = blockIdx.x * blockDim.x + threadIdx.x;
+  if (row_block_id > num_row_blocks) {
+    // No synchronization is performed in this kernel, thus it is safe to return
+    return;
+  }
+  if (row_block_id == num_row_blocks) {
+    // one extra thread sets the first element
+    rows[0] = 0;
+    return;
+  }
+  const auto& row_block = row_blocks[row_block_id];
+  int row_nnz = 0;
+  const auto first_cell = cells + first_cell_in_row_block[row_block_id];
+  const auto last_cell = cells + first_cell_in_row_block[row_block_id + 1];
+  for (auto cell = first_cell; cell < last_cell; ++cell) {
+    row_nnz += col_blocks[cell->block_id].size;
+  }
+  const int first_row = row_block.position;
+  const int last_row = first_row + row_block.size;
+  for (int i = first_row; i < last_row; ++i) {
+    rows[i + 1] = row_nnz;
+    row_block_ids[i] = row_block_id;
+  }
+}
+
+// Row-wise creation of CRS structure
+// Inputs:
+// - num_rows: number of rows in matrix
+// - first_cell_in_row_block: index of the first cell of the row-block; size:
+// num_row_blocks + 1
+// - cells: cells of block-sparse structure as a continuous array
+// - row_blocks: row blocks of block-sparse structure stored sequentially
+// - col_blocks: column blocks of block-sparse structure stored sequentially
+// - row_block_ids: index of row-block that corresponds to row
+// - rows: row-index array of CRS structure
+// Outputs:
+// - cols: column-index array of CRS structure
+// Computaion is perform row-wise
+__global__ void ComputeColumnsAndPermutation(
+    int num_rows,
+    const int* __restrict__ first_cell_in_row_block,
+    const Cell* __restrict__ cells,
+    const Block* __restrict__ row_blocks,
+    const Block* __restrict__ col_blocks,
+    const int* __restrict__ row_block_ids,
+    const int* __restrict__ rows,
+    int* __restrict__ cols) {
+  const int row = blockIdx.x * blockDim.x + threadIdx.x;
+  if (row >= num_rows) {
+    // No synchronization is performed in this kernel, thus it is safe to return
+    return;
+  }
+  const int row_block_id = row_block_ids[row];
+  // position in crs matrix
+  int crs_position = rows[row];
+  const auto first_cell = cells + first_cell_in_row_block[row_block_id];
+  const auto last_cell = cells + first_cell_in_row_block[row_block_id + 1];
+  // For reach cell of row-block only current row is being filled
+  for (auto cell = first_cell; cell < last_cell; ++cell) {
+    const auto& col_block = col_blocks[cell->block_id];
+    const int col_block_size = col_block.size;
+    int column_idx = col_block.position;
+    // Column indices for each element of row_in_block row of current cell
+    for (int i = 0; i < col_block_size; ++i, ++crs_position) {
+      cols[crs_position] = column_idx++;
+    }
+  }
+}
+
+void FillCRSStructure(const int num_row_blocks,
+                      const int num_rows,
+                      const int* first_cell_in_row_block,
+                      const Cell* cells,
+                      const Block* row_blocks,
+                      const Block* col_blocks,
+                      int* rows,
+                      int* cols,
+                      cudaStream_t stream) {
+  // Set number of non-zeros per row in rows array and row to row-block map in
+  // row_block_ids array
+  int* row_block_ids;
+  cudaMallocAsync(&row_block_ids, sizeof(int) * num_rows, stream);
+  const int num_blocks_blockwise = NumBlocksInGrid(num_row_blocks + 1);
+  RowBlockIdAndNNZ<<<num_blocks_blockwise, kCudaBlockSize, 0, stream>>>(
+      num_row_blocks,
+      first_cell_in_row_block,
+      cells,
+      row_blocks,
+      col_blocks,
+      rows,
+      row_block_ids);
+  // Finalize row-index array of CRS strucure by computing prefix sum
+  thrust::inclusive_scan(
+      thrust::cuda::par_nosync.on(stream), rows, rows + num_rows + 1, rows);
+
+  // Fill cols array of CRS structure and permutation from block-sparse to CRS
+  const int num_blocks_rowwise = NumBlocksInGrid(num_rows);
+  ComputeColumnsAndPermutation<<<num_blocks_rowwise,
+                                 kCudaBlockSize,
+                                 0,
+                                 stream>>>(num_rows,
+                                           first_cell_in_row_block,
+                                           cells,
+                                           row_blocks,
+                                           col_blocks,
+                                           row_block_ids,
+                                           rows,
+                                           cols);
+  cudaFreeAsync(row_block_ids, stream);
+}
+
+template <typename T, typename Predicate>
+__device__ int PartitionPoint(const T* data,
+                              int first,
+                              int last,
+                              Predicate&& predicate) {
+  if (!predicate(data[first])) {
+    return first;
+  }
+  while (last - first > 1) {
+    const auto midpoint = first + (last - first) / 2;
+    if (predicate(data[midpoint])) {
+      first = midpoint;
+    } else {
+      last = midpoint;
+    }
+  }
+  return last;
+}
+
+// Element-wise reordering of block-sparse values
+// - first_cell_pos - position of the first cell of corresponding sub-matrix;
+// only used if use_first_cell is true; otherwise the first cell in row-block is
+// used
+// - block_sparse_values - segment of block-sparse values starting from
+// block_sparse_offset, containing num_values
+__global__ void PermuteToCrsKernel(
+    const int block_sparse_offset,
+    const int num_values,
+    const int num_row_blocks,
+    const int* __restrict__ first_cell_in_row_block,
+    const Cell* __restrict__ cells,
+    const Block* __restrict__ row_blocks,
+    const Block* __restrict__ col_blocks,
+    const int* __restrict__ crs_rows,
+    const double* __restrict__ block_sparse_values,
+    double* __restrict__ crs_values) {
+  const int value_id = blockIdx.x * blockDim.x + threadIdx.x;
+  if (value_id >= num_values) {
+    return;
+  }
+  const int block_sparse_value_id = value_id + block_sparse_offset;
+  // Find the corresponding row-block with a binary search
+  const int row_block_id =
+      PartitionPoint(
+          first_cell_in_row_block,
+          0,
+          num_row_blocks,
+          [cells,
+           block_sparse_value_id] __device__(const int row_block_offset) {
+            return cells[row_block_offset].position <= block_sparse_value_id;
+          }) -
+      1;
+  // Find cell and calculate offset within the row with a linear scan
+  const auto& row_block = row_blocks[row_block_id];
+  const auto first_cell = cells + first_cell_in_row_block[row_block_id];
+  const auto last_cell = cells + first_cell_in_row_block[row_block_id + 1];
+  const int row_block_size = row_block.size;
+  int num_cols_before = 0;
+  for (const Cell* cell = first_cell; cell < last_cell; ++cell) {
+    const auto& col_block = col_blocks[cell->block_id];
+    const int col_block_size = col_block.size;
+    const int cell_size = row_block_size * col_block_size;
+    if (cell->position + cell_size > block_sparse_value_id) {
+      const int pos_in_cell = block_sparse_value_id - cell->position;
+      const int row_in_cell = pos_in_cell / col_block_size;
+      const int col_in_cell = pos_in_cell % col_block_size;
+      const int row = row_in_cell + row_block.position;
+      crs_values[crs_rows[row] + num_cols_before + col_in_cell] =
+          block_sparse_values[value_id];
+      break;
+    }
+    num_cols_before += col_block_size;
+  }
+}
+
+void PermuteToCRS(const int block_sparse_offset,
+                  const int num_values,
+                  const int num_row_blocks,
+                  const int* first_cell_in_row_block,
+                  const Cell* cells,
+                  const Block* row_blocks,
+                  const Block* col_blocks,
+                  const int* crs_rows,
+                  const double* block_sparse_values,
+                  double* crs_values,
+                  cudaStream_t stream) {
+  const int num_blocks_valuewise = NumBlocksInGrid(num_values);
+  PermuteToCrsKernel<<<num_blocks_valuewise, kCudaBlockSize, 0, stream>>>(
+      block_sparse_offset,
+      num_values,
+      num_row_blocks,
+      first_cell_in_row_block,
+      cells,
+      row_blocks,
+      col_blocks,
+      crs_rows,
+      block_sparse_values,
+      crs_values);
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/cuda_kernels_bsm_to_crs.h b/internal/ceres/cuda_kernels_bsm_to_crs.h
new file mode 100644
index 0000000..b6da69b
--- /dev/null
+++ b/internal/ceres/cuda_kernels_bsm_to_crs.h
@@ -0,0 +1,76 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 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.
+//
+// Authors: dmitriy.korchemkin@gmail.com (Dmitriy Korchemkin)
+
+#ifndef CERES_INTERNAL_CUDA_KERNELS_BSM_TO_CRS_H_
+#define CERES_INTERNAL_CUDA_KERNELS_BSM_TO_CRS_H_
+
+#include "ceres/internal/config.h"
+
+#ifndef CERES_NO_CUDA
+
+#include "cuda_runtime.h"
+
+namespace ceres {
+namespace internal {
+class Block;
+class Cell;
+
+// Compute structure of CRS matrix using block-sparse structure.
+// Arrays corresponding to CRS matrix are to be allocated by caller
+void FillCRSStructure(const int num_row_blocks,
+                      const int num_rows,
+                      const int* first_cell_in_row_block,
+                      const Cell* cells,
+                      const Block* row_blocks,
+                      const Block* col_blocks,
+                      int* rows,
+                      int* cols,
+                      cudaStream_t stream);
+
+// Permute segment of values from block-sparse matrix with sequential layout to
+// CRS order. Segment starts at block_sparse_offset and has length of num_values
+void PermuteToCRS(const int block_sparse_offset,
+                  const int num_values,
+                  const int num_row_blocks,
+                  const int* first_cell_in_row_block,
+                  const Cell* cells,
+                  const Block* row_blocks,
+                  const Block* col_blocks,
+                  const int* crs_rows,
+                  const double* block_sparse_values,
+                  double* crs_values,
+                  cudaStream_t stream);
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_NO_CUDA
+
+#endif  // CERES_INTERNAL_CUDA_KERNELS_BSM_TO_CRS_H_
diff --git a/internal/ceres/cuda_kernels_utils.h b/internal/ceres/cuda_kernels_utils.h
new file mode 100644
index 0000000..45d5d3b
--- /dev/null
+++ b/internal/ceres/cuda_kernels_utils.h
@@ -0,0 +1,56 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: joydeepb@cs.utexas.edu (Joydeep Biswas)
+
+#ifndef CERES_INTERNAL_CUDA_KERNELS_UTILS_H_
+#define CERES_INTERNAL_CUDA_KERNELS_UTILS_H_
+
+namespace ceres {
+namespace internal {
+
+// Parallel execution on CUDA device requires splitting job into blocks of a
+// fixed size. We use block-size of kCudaBlockSize for all kernels that do not
+// require any specific block size. 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;
+
+// Compute number of blocks of kCudaBlockSize that span over 1-d grid with
+// dimension size. Note that 1-d grid dimension is limited by 2^31-1 in CUDA,
+// thus a signed int is used as an argument.
+inline int NumBlocksInGrid(int size) {
+  return (size + kCudaBlockSize - 1) / kCudaBlockSize;
+}
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_CUDA_KERNELS_UTILS_H_
diff --git a/internal/ceres/cuda_kernels_vector_ops.cu.cc b/internal/ceres/cuda_kernels_vector_ops.cu.cc
new file mode 100644
index 0000000..b5a8e40
--- /dev/null
+++ b/internal/ceres/cuda_kernels_vector_ops.cu.cc
@@ -0,0 +1,123 @@
+// 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/cuda_kernels_vector_ops.h"
+
+#include <cuda_runtime.h>
+
+#include "ceres/cuda_kernels_utils.h"
+
+namespace ceres {
+namespace internal {
+
+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 = NumBlocksInGrid(size);
+  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 = NumBlocksInGrid(size);
+  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 = NumBlocksInGrid(size);
+  SetZeroKernel<float><<<num_blocks, kCudaBlockSize, 0, stream>>>(output, size);
+}
+
+void CudaSetZeroFP64(double* output, const int size, cudaStream_t stream) {
+  const int num_blocks = NumBlocksInGrid(size);
+  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 = NumBlocksInGrid(size);
+  XPlusEqualsYKernel<float, double>
+      <<<num_blocks, kCudaBlockSize, 0, stream>>>(x, y, size);
+}
+
+__global__ void CudaDtDxpyKernel(double* __restrict__ y,
+                                 const double* D,
+                                 const double* __restrict__ x,
+                                 const int size) {
+  const int i = blockIdx.x * blockDim.x + threadIdx.x;
+  if (i < size) {
+    y[i] = y[i] + D[i] * D[i] * x[i];
+  }
+}
+
+void CudaDtDxpy(double* y,
+                const double* D,
+                const double* x,
+                const int size,
+                cudaStream_t stream) {
+  const int num_blocks = NumBlocksInGrid(size);
+  CudaDtDxpyKernel<<<num_blocks, kCudaBlockSize, 0, stream>>>(y, D, x, size);
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/cuda_kernels.h b/internal/ceres/cuda_kernels_vector_ops.h
similarity index 66%
rename from internal/ceres/cuda_kernels.h
rename to internal/ceres/cuda_kernels_vector_ops.h
index 61f945a..b179201 100644
--- a/internal/ceres/cuda_kernels.h
+++ b/internal/ceres/cuda_kernels_vector_ops.h
@@ -28,8 +28,8 @@
 //
 // Author: joydeepb@cs.utexas.edu (Joydeep Biswas)
 
-#ifndef CERES_INTERNAL_CUDA_KERNELS_H_
-#define CERES_INTERNAL_CUDA_KERNELS_H_
+#ifndef CERES_INTERNAL_CUDA_KERNELS_VECTOR_OPS_H_
+#define CERES_INTERNAL_CUDA_KERNELS_VECTOR_OPS_H_
 
 #include "ceres/internal/config.h"
 
@@ -75,39 +75,9 @@
                 const int size,
                 cudaStream_t stream);
 
-// Compute structure of CRS matrix and permutation of values using block-sparse
-// structure and temporary array row_block_ids. Array row_block_ids of size
-// num_rows will be filled with indices of row-blocks corresponding to rows of
-// CRS matrix. Arrays corresponding to CRS matrix, permutation and row_block_ids
-// arrays are to be allocated by caller
-void FillCRSStructure(const int num_row_blocks,
-                      const int num_rows,
-                      const int* row_block_offsets,
-                      const Cell* cells,
-                      const Block* row_blocks,
-                      const Block* col_blocks,
-                      int* rows,
-                      int* cols,
-                      int* row_block_ids,
-                      int* permutation,
-                      cudaStream_t stream);
-
-// Permute block of block-sparse values using permutation
-// Pointer block_sparse_values corresponds to a block of num_values values from
-// block-sparse matrix at the offset from begining. Pointer output corresponds
-// to values of CRS matrix. Array permutation stores permutation from
-// block-sparse to CRS matrix with permutation[i] being an index of i-th value
-// of block-sparse matrix in values of CRS matrix
-void PermuteValues(const int offset,
-                   const int num_values,
-                   const int* permutation,
-                   const double* block_sparse_values,
-                   double* crs_values,
-                   cudaStream_t stream);
-
 }  // namespace internal
 }  // namespace ceres
 
 #endif  // CERES_NO_CUDA
 
-#endif  // CERES_INTERNAL_CUDA_KERNELS_H_
+#endif  // CERES_INTERNAL_CUDA_KERNELS_VECTOR_OPS_H_
diff --git a/internal/ceres/cuda_kernels_test.cc b/internal/ceres/cuda_kernels_vector_ops_test.cc
similarity index 99%
rename from internal/ceres/cuda_kernels_test.cc
rename to internal/ceres/cuda_kernels_vector_ops_test.cc
index 41036a3..7f6e229 100644
--- a/internal/ceres/cuda_kernels_test.cc
+++ b/internal/ceres/cuda_kernels_vector_ops_test.cc
@@ -28,7 +28,7 @@
 //
 // Author: joydeepb@cs.utexas.edu (Joydeep Biswas)
 
-#include "ceres/cuda_kernels.h"
+#include "ceres/cuda_kernels_vector_ops.h"
 
 #include <math.h>
 
diff --git a/internal/ceres/cuda_sparse_matrix.cc b/internal/ceres/cuda_sparse_matrix.cc
index 0ae8c5d..e9ca9cc 100644
--- a/internal/ceres/cuda_sparse_matrix.cc
+++ b/internal/ceres/cuda_sparse_matrix.cc
@@ -52,7 +52,7 @@
 #ifndef CERES_NO_CUDA
 
 #include "ceres/cuda_buffer.h"
-#include "ceres/cuda_kernels.h"
+#include "ceres/cuda_kernels_vector_ops.h"
 #include "ceres/cuda_vector.h"
 #include "cuda_runtime_api.h"
 #include "cusparse.h"
diff --git a/internal/ceres/cuda_sparse_matrix.h b/internal/ceres/cuda_sparse_matrix.h
index ad55170..f5fcb91 100644
--- a/internal/ceres/cuda_sparse_matrix.h
+++ b/internal/ceres/cuda_sparse_matrix.h
@@ -82,6 +82,10 @@
   int num_cols() const { return num_cols_; }
   int num_nonzeros() const { return num_nonzeros_; }
 
+  const int32_t* rows() const { return rows_.data(); }
+  const int32_t* cols() const { return cols_.data(); }
+  const double* values() const { return values_.data(); }
+
   int32_t* mutable_rows() { return rows_.data(); }
   int32_t* mutable_cols() { return cols_.data(); }
   double* mutable_values() { return values_.data(); }
diff --git a/internal/ceres/cuda_vector.cc b/internal/ceres/cuda_vector.cc
index e2143de..d434f36 100644
--- a/internal/ceres/cuda_vector.cc
+++ b/internal/ceres/cuda_vector.cc
@@ -44,7 +44,7 @@
 #ifndef CERES_NO_CUDA
 
 #include "ceres/cuda_buffer.h"
-#include "ceres/cuda_kernels.h"
+#include "ceres/cuda_kernels_vector_ops.h"
 #include "ceres/cuda_vector.h"
 #include "cublas_v2.h"
 
diff --git a/internal/ceres/cuda_vector.h b/internal/ceres/cuda_vector.h
index 34f3947..46661cf 100644
--- a/internal/ceres/cuda_vector.h
+++ b/internal/ceres/cuda_vector.h
@@ -50,7 +50,7 @@
 #ifndef CERES_NO_CUDA
 
 #include "ceres/cuda_buffer.h"
-#include "ceres/cuda_kernels.h"
+#include "ceres/cuda_kernels_vector_ops.h"
 #include "ceres/internal/eigen.h"
 #include "cublas_v2.h"
 #include "cusparse.h"
diff --git a/internal/ceres/dense_cholesky.cc b/internal/ceres/dense_cholesky.cc
index f19ca7e..fb95536 100644
--- a/internal/ceres/dense_cholesky.cc
+++ b/internal/ceres/dense_cholesky.cc
@@ -41,7 +41,7 @@
 
 #ifndef CERES_NO_CUDA
 #include "ceres/context_impl.h"
-#include "ceres/cuda_kernels.h"
+#include "ceres/cuda_kernels_vector_ops.h"
 #include "cuda_runtime.h"
 #include "cusolverDn.h"
 #endif  // CERES_NO_CUDA
diff --git a/internal/ceres/evaluation_benchmark.cc b/internal/ceres/evaluation_benchmark.cc
index 2ba03ed..c23c881 100644
--- a/internal/ceres/evaluation_benchmark.cc
+++ b/internal/ceres/evaluation_benchmark.cc
@@ -96,7 +96,7 @@
   }
 
   std::unique_ptr<BlockSparseMatrix> CreateBlockSparseJacobian(
-      ContextImpl* context) {
+      ContextImpl* context, bool sequential) {
     auto problem = bal_problem->mutable_problem();
     auto problem_impl = problem->mutable_impl();
     CHECK(problem_impl != nullptr);
@@ -116,6 +116,30 @@
     auto block_sparse = downcast_unique_ptr<BlockSparseMatrix>(jacobian);
     CHECK(block_sparse != nullptr);
 
+    if (sequential) {
+      auto block_structure_sequential =
+          std::make_unique<CompressedRowBlockStructure>(
+              *block_sparse->block_structure());
+      int num_nonzeros = 0;
+      for (auto& row_block : block_structure_sequential->rows) {
+        const int row_block_size = row_block.block.size;
+        for (auto& cell : row_block.cells) {
+          const int col_block_size =
+              block_structure_sequential->cols[cell.block_id].size;
+          cell.position = num_nonzeros;
+          num_nonzeros += col_block_size * row_block_size;
+        }
+      }
+      block_sparse = std::make_unique<BlockSparseMatrix>(
+          block_structure_sequential.release(),
+#ifndef CERES_NO_CUDA
+          true
+#else
+          false
+#endif
+      );
+    }
+
     std::mt19937 rng;
     std::normal_distribution<double> rnorm;
     const int nnz = block_sparse->num_nonzeros();
@@ -123,6 +147,7 @@
     for (int i = 0; i < nnz; ++i) {
       values[i] = rnorm(rng);
     }
+
     return block_sparse;
   }
 
@@ -134,11 +159,20 @@
 
   const BlockSparseMatrix* BlockSparseJacobian(ContextImpl* context) {
     if (!block_sparse_jacobian) {
-      block_sparse_jacobian = CreateBlockSparseJacobian(context);
+      block_sparse_jacobian = CreateBlockSparseJacobian(context, true);
     }
     return block_sparse_jacobian.get();
   }
 
+  const BlockSparseMatrix* BlockSparseJacobianPartitioned(
+      ContextImpl* context) {
+    if (!block_sparse_jacobian_partitioned) {
+      block_sparse_jacobian_partitioned =
+          CreateBlockSparseJacobian(context, false);
+    }
+    return block_sparse_jacobian_partitioned.get();
+  }
+
   const CompressedRowSparseMatrix* CompressedRowSparseJacobian(
       ContextImpl* context) {
     if (!crs_jacobian) {
@@ -149,7 +183,7 @@
 
   std::unique_ptr<PartitionedView> PartitionedMatrixViewJacobian(
       const LinearSolver::Options& options) {
-    auto block_sparse = BlockSparseJacobian(options.context);
+    auto block_sparse = BlockSparseJacobianPartitioned(options.context);
     return std::make_unique<PartitionedView>(options, *block_sparse);
   }
 
@@ -171,7 +205,7 @@
 
   const ImplicitSchurComplement* ImplicitSchurComplementWithoutDiagonal(
       const LinearSolver::Options& options) {
-    auto block_sparse = BlockSparseJacobian(options.context);
+    auto block_sparse = BlockSparseJacobianPartitioned(options.context);
     implicit_schur_complement =
         std::make_unique<ImplicitSchurComplement>(options);
     implicit_schur_complement->Init(*block_sparse, nullptr, b.data());
@@ -180,7 +214,7 @@
 
   const ImplicitSchurComplement* ImplicitSchurComplementWithDiagonal(
       const LinearSolver::Options& options) {
-    auto block_sparse = BlockSparseJacobian(options.context);
+    auto block_sparse = BlockSparseJacobianPartitioned(options.context);
     implicit_schur_complement_diag =
         std::make_unique<ImplicitSchurComplement>(options);
     implicit_schur_complement_diag->Init(*block_sparse, D.data(), b.data());
@@ -192,6 +226,7 @@
   Vector b;
   std::unique_ptr<BundleAdjustmentProblem> bal_problem;
   std::unique_ptr<PreprocessedProblem> preprocessed_problem;
+  std::unique_ptr<BlockSparseMatrix> block_sparse_jacobian_partitioned;
   std::unique_ptr<BlockSparseMatrix> block_sparse_jacobian;
   std::unique_ptr<CompressedRowSparseMatrix> crs_jacobian;
   std::unique_ptr<BlockSparseMatrix> block_diagonal_ete;