diff --git a/BUILD b/BUILD
index a82d2cf..504a23e 100644
--- a/BUILD
+++ b/BUILD
@@ -141,6 +141,7 @@
     "solver",
     "sparse_cholesky",
     "sparse_normal_cholesky_solver",
+    "subset_preconditioner",
     "system",
     "tiny_solver_autodiff_function",
     "tiny_solver_cost_function_adapter",
diff --git a/bazel/ceres.bzl b/bazel/ceres.bzl
index 833c937..d9d7502 100644
--- a/bazel/ceres.bzl
+++ b/bazel/ceres.bzl
@@ -111,6 +111,7 @@
     "sparse_normal_cholesky_solver.cc",
     "split.cc",
     "stringprintf.cc",
+    "subset_preconditioner.cc",
     "suitesparse.cc",
     "thread_token_provider.cc",
     "triplet_sparse_matrix.cc",
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 884fd57..91ffc11 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -111,6 +111,7 @@
     sparse_matrix.cc
     sparse_cholesky.cc
     sparse_normal_cholesky_solver.cc
+    subset_preconditioner.cc
     split.cc
     stringprintf.cc
     suitesparse.cc
@@ -359,6 +360,7 @@
   ceres_test(solver)
   ceres_test(sparse_cholesky)
   ceres_test(sparse_normal_cholesky_solver)
+  ceres_test(subset_preconditioner)
   ceres_test(system)
   ceres_test(tiny_solver)
   ceres_test(tiny_solver_autodiff_function)
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index 214daef..c814a72 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -280,9 +280,12 @@
 }
 
 void BlockSparseMatrix::AppendRows(const BlockSparseMatrix& m) {
+  CHECK_EQ(m.num_cols(), num_cols());
+  const CompressedRowBlockStructure* m_bs = m.block_structure();
+  CHECK_EQ(m_bs->cols.size(), block_structure_->cols.size());
+
   const int old_num_nonzeros = num_nonzeros_;
   const int old_num_row_blocks = block_structure_->rows.size();
-  const CompressedRowBlockStructure* m_bs = m.block_structure();
   block_structure_->rows.resize(old_num_row_blocks + m_bs->rows.size());
 
   for (int i = 0; i < m_bs->rows.size(); ++i) {
@@ -336,29 +339,33 @@
   CHECK_GT(options.min_row_block_size, 0);
   CHECK_GT(options.max_row_block_size, 0);
   CHECK_LE(options.min_row_block_size, options.max_row_block_size);
-  CHECK_GT(options.num_col_blocks, 0);
-  CHECK_GT(options.min_col_block_size, 0);
-  CHECK_GT(options.max_col_block_size, 0);
-  CHECK_LE(options.min_col_block_size, options.max_col_block_size);
   CHECK_GT(options.block_density, 0.0);
   CHECK_LE(options.block_density, 1.0);
 
   CompressedRowBlockStructure* bs = new CompressedRowBlockStructure();
-    // Generate the col block structure.
-  int col_block_position = 0;
-  for (int i = 0; i < options.num_col_blocks; ++i) {
-    // Generate a random integer in [min_col_block_size, max_col_block_size]
-    const int delta_block_size =
-        Uniform(options.max_col_block_size - options.min_col_block_size);
-    const int col_block_size = options.min_col_block_size + delta_block_size;
-    bs->cols.push_back(Block(col_block_size, col_block_position));
-    col_block_position += col_block_size;
-  }
+  if (options.col_blocks.empty()) {
+    CHECK_GT(options.num_col_blocks, 0);
+    CHECK_GT(options.min_col_block_size, 0);
+    CHECK_GT(options.max_col_block_size, 0);
+    CHECK_LE(options.min_col_block_size, options.max_col_block_size);
 
+    // Generate the col block structure.
+    int col_block_position = 0;
+    for (int i = 0; i < options.num_col_blocks; ++i) {
+      // Generate a random integer in [min_col_block_size, max_col_block_size]
+      const int delta_block_size =
+          Uniform(options.max_col_block_size - options.min_col_block_size);
+      const int col_block_size = options.min_col_block_size + delta_block_size;
+      bs->cols.push_back(Block(col_block_size, col_block_position));
+      col_block_position += col_block_size;
+    }
+  } else {
+    bs->cols = options.col_blocks;
+  }
 
   bool matrix_has_blocks = false;
   while (!matrix_has_blocks) {
-    LOG(INFO) << "clearing";
+    VLOG(1) << "Clearing";
     bs->rows.clear();
     int row_block_position = 0;
     int value_position = 0;
@@ -372,7 +379,7 @@
       row.block.size = row_block_size;
       row.block.position = row_block_position;
       row_block_position += row_block_size;
-      for (int c = 0; c < options.num_col_blocks; ++c) {
+      for (int c = 0; c < bs->cols.size(); ++c) {
         if (RandDouble() > options.block_density) continue;
 
         row.cells.push_back(Cell());
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index 17491bf..abde2a6 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -117,6 +117,11 @@
     // present in the matrix. A given random matrix will not have
     // precisely this density.
     double block_density;
+
+    // If col_blocks is non-empty, then the generated random matrix
+    // has this block structure and the column related options in this
+    // struct are ignored.
+    std::vector<Block> col_blocks;
   };
 
   // Create a random BlockSparseMatrix whose entries are normally
diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h
index a248eae..5d293b8 100644
--- a/internal/ceres/preconditioner.h
+++ b/internal/ceres/preconditioner.h
@@ -51,6 +51,8 @@
         : type(JACOBI),
           visibility_clustering_type(CANONICAL_VIEWS),
           sparse_linear_algebra_library_type(SUITE_SPARSE),
+          use_postordering(false),
+          subset_preconditioner_start_row_block(-1),
           num_threads(1),
           row_block_size(Eigen::Dynamic),
           e_block_size(Eigen::Dynamic),
@@ -61,6 +63,21 @@
     VisibilityClusteringType visibility_clustering_type;
     SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
 
+    // When using the subset preconditioner, all row blocks starting
+    // from this row block are used to construct the preconditioner.
+    //
+    // i.e., the Jacobian matrix A is horizonatally partitioned as
+    //
+    // A = [P]
+    //     [Q]
+    //
+    // where P has subset_preconditioner_start_row_block row blocks,
+    // and the preconditioner is the inverse of the matrix Q'Q.
+    int subset_preconditioner_start_row_block;
+
+    // See solver.h for information about these flags.
+    bool use_postordering;
+
     // If possible, how many threads the preconditioner can use.
     int num_threads;
 
diff --git a/internal/ceres/subset_preconditioner.cc b/internal/ceres/subset_preconditioner.cc
new file mode 100644
index 0000000..e970b91
--- /dev/null
+++ b/internal/ceres/subset_preconditioner.cc
@@ -0,0 +1,112 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 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: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/subset_preconditioner.h"
+
+#include <string>
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/inner_product_computer.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/linear_solver.h"
+#include "ceres/sparse_cholesky.h"
+#include "ceres/types.h"
+
+namespace ceres {
+namespace internal {
+
+SubsetPreconditioner::SubsetPreconditioner(
+    const Preconditioner::Options& options, const BlockSparseMatrix& A)
+    : options_(options), num_cols_(A.num_cols()) {
+  sparse_cholesky_.reset(
+      SparseCholesky::Create(options_.sparse_linear_algebra_library_type,
+                             options_.use_postordering ? AMD : NATURAL));
+  CHECK_GE(options_.subset_preconditioner_start_row_block, 0);
+}
+
+SubsetPreconditioner::~SubsetPreconditioner() {}
+
+void SubsetPreconditioner::RightMultiply(const double* x, double* y) const {
+  CHECK_NOTNULL(x);
+  CHECK_NOTNULL(y);
+  std::string message;
+  sparse_cholesky_->Solve(x, y, &message);
+}
+
+bool SubsetPreconditioner::UpdateImpl(const BlockSparseMatrix& A,
+                                      const double* D) {
+  BlockSparseMatrix* m = const_cast<BlockSparseMatrix*>(&A);
+  const CompressedRowBlockStructure* bs = m->block_structure();
+
+  // A = [P]
+  //     [Q]
+
+  // Now add D to A if needed.
+  if (D != NULL) {
+    // A = [P]
+    //     [Q]
+    //     [D]
+    scoped_ptr<BlockSparseMatrix> regularizer(
+        BlockSparseMatrix::CreateDiagonalMatrix(D, bs->cols));
+    m->AppendRows(*regularizer);
+  }
+
+  if (inner_product_computer_.get() == NULL) {
+    inner_product_computer_.reset(InnerProductComputer::Create(
+        *m,
+        options_.subset_preconditioner_start_row_block,
+        bs->rows.size(),
+        sparse_cholesky_->StorageType()));
+  }
+
+  // Compute inner_product = [Q'*Q + D'*D]
+  inner_product_computer_->Compute();
+
+  // Unappend D if needed.
+  if (D != NULL) {
+    // A = [P]
+    //     [Q]
+    m->DeleteRowBlocks(bs->cols.size());
+  }
+
+  std::string message;
+  // Compute L. s.t., LL' = Q'*Q + D'*D
+  const LinearSolverTerminationType termination_type =
+      sparse_cholesky_->Factorize(inner_product_computer_->mutable_result(),
+                                  &message);
+  if (termination_type != LINEAR_SOLVER_SUCCESS) {
+    LOG(ERROR) << "Preconditioner factorization failed: " << message;
+    return false;
+  }
+
+  return true;
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/subset_preconditioner.h b/internal/ceres/subset_preconditioner.h
new file mode 100644
index 0000000..062253b
--- /dev/null
+++ b/internal/ceres/subset_preconditioner.h
@@ -0,0 +1,91 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 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: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_SUBSET_PRECONDITIONER_H_
+#define CERES_INTERNAL_SUBSET_PRECONDITIONER_H_
+
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/preconditioner.h"
+
+namespace ceres {
+namespace internal {
+
+class BlockSparseMatrix;
+class SparseCholesky;
+class InnerProductComputer;
+
+// Subset preconditioning, uses a subset of the rows of the Jacobian
+// to construct a preconditioner for the normal equations.
+//
+// To keep the interface simple, we assume that the matrix A has
+// already been re-ordered that the user wishes to some subset of the
+// bottom row blocks of the matrix as the preconditioner. This is
+// controlled by
+// Preconditioner::Options::subset_preconditioner_start_row_block.
+//
+// When using the subset preconditioner, all row blocks starting
+// from this row block are used to construct the preconditioner.
+//
+// More precisely the matrix A is horizontally partitioned as
+//
+// A = [P]
+//     [Q]
+//
+// where P as subset_preconditioner_start_row_block row blocks,
+// and the preconditioner is the inverse of the matrix Q'Q.
+//
+// Obviously, the smaller this number, the more accurate and
+// computationally expensive this preconditioner will be.
+//
+// See the tests for example usage.
+class SubsetPreconditioner : public BlockSparseMatrixPreconditioner {
+ public:
+  SubsetPreconditioner(const Preconditioner::Options& options,
+                       const BlockSparseMatrix& A);
+  virtual ~SubsetPreconditioner();
+
+  // Preconditioner interface
+  virtual void RightMultiply(const double* x, double* y) const;
+  virtual int num_rows() const { return num_cols_; }
+  virtual int num_cols() const { return num_cols_; }
+
+ private:
+  virtual bool UpdateImpl(const BlockSparseMatrix& A, const double* D);
+
+  const Preconditioner::Options options_;
+  const int num_cols_;
+  scoped_ptr<SparseCholesky> sparse_cholesky_;
+  scoped_ptr<InnerProductComputer> inner_product_computer_;
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_SUBSET_PRECONDITIONER_H_
diff --git a/internal/ceres/subset_preconditioner_test.cc b/internal/ceres/subset_preconditioner_test.cc
new file mode 100644
index 0000000..63dbbaf
--- /dev/null
+++ b/internal/ceres/subset_preconditioner_test.cc
@@ -0,0 +1,194 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 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: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/subset_preconditioner.h"
+#include "Eigen/Dense"
+#include "Eigen/SparseCore"
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/inner_product_computer.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+// TODO(sameeragarwal): Refactor the following two functions out of
+// here and sparse_cholesky_test.cc into a more suitable place.
+template <Eigen::UpLoType UpLoType>
+bool SolveLinearSystemUsingEigen(const Matrix& lhs,
+                                 const Vector rhs,
+                                 Vector* solution) {
+  Eigen::LLT<Matrix, UpLoType> llt = lhs.selfadjointView<UpLoType>().llt();
+  if (llt.info() != Eigen::Success) {
+    return false;
+  }
+  *solution = llt.solve(rhs);
+  return (llt.info() == Eigen::Success);
+}
+
+// Use Eigen's Dense Cholesky solver to compute the solution to a
+// sparse linear system.
+bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
+                             const Vector& rhs,
+                             Vector* solution) {
+  Matrix dense_triangular_lhs;
+  lhs.ToDenseMatrix(&dense_triangular_lhs);
+  if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+    Matrix full_lhs = dense_triangular_lhs.selfadjointView<Eigen::Upper>();
+    return SolveLinearSystemUsingEigen<Eigen::Upper>(full_lhs, rhs, solution);
+  }
+  return SolveLinearSystemUsingEigen<Eigen::Lower>(
+      dense_triangular_lhs, rhs, solution);
+}
+
+typedef ::testing::tuple<SparseLinearAlgebraLibraryType, bool> Param;
+
+std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+  Param param = info.param;
+  std::stringstream ss;
+  ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
+     << (::testing::get<1>(param) ? "Diagonal" : "NoDiagonal");
+  return ss.str();
+}
+
+class SubsetPreconditionerTest : public ::testing::TestWithParam<Param> {
+ protected:
+  virtual void SetUp() {
+    BlockSparseMatrix::RandomMatrixOptions options;
+    options.num_col_blocks = 4;
+    options.min_col_block_size = 1;
+    options.max_col_block_size = 4;
+    options.num_row_blocks = 8;
+    options.min_row_block_size = 1;
+    options.max_row_block_size = 4;
+    options.block_density = 0.9;
+
+    m_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
+    start_row_block_ = m_->block_structure()->rows.size();
+
+    // Ensure that the bottom part of the matrix has the same column
+    // block structure.
+    options.col_blocks = m_->block_structure()->cols;
+    b_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
+    m_->AppendRows(*b_);
+
+    // Create a Identity block diagonal matrix with the same column
+    // block structure.
+    diagonal_ = Vector::Ones(m_->num_cols());
+    block_diagonal_.reset(BlockSparseMatrix::CreateDiagonalMatrix(
+        diagonal_.data(), b_->block_structure()->cols));
+
+    // Unconditionally add the block diagonal to the matrix b_,
+    // because either it is either part of b_ to make it full rank, or
+    // we pass the same diagonal matrix later as the parameter D. In
+    // either case the preconditioner matrix is b_' b + D'D.
+    b_->AppendRows(*block_diagonal_);
+    inner_product_computer_.reset(InnerProductComputer::Create(
+        *b_, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+    inner_product_computer_->Compute();
+  }
+
+  scoped_ptr<BlockSparseMatrix> m_;
+  scoped_ptr<BlockSparseMatrix> b_;
+  scoped_ptr<BlockSparseMatrix> block_diagonal_;
+  scoped_ptr<InnerProductComputer> inner_product_computer_;
+  scoped_ptr<Preconditioner> preconditioner_;
+  Vector diagonal_;
+  int start_row_block_;
+};
+
+TEST_P(SubsetPreconditionerTest, foo) {
+  Param param = GetParam();
+  Preconditioner::Options options;
+  options.subset_preconditioner_start_row_block = start_row_block_;
+  options.sparse_linear_algebra_library_type = ::testing::get<0>(param);
+  preconditioner_.reset(new SubsetPreconditioner(options, *m_));
+
+  const bool with_diagonal = ::testing::get<1>(param);
+  if (!with_diagonal) {
+    m_->AppendRows(*block_diagonal_);
+  }
+
+  EXPECT_TRUE(
+      preconditioner_->Update(*m_, with_diagonal ? diagonal_.data() : NULL));
+
+  // Repeatedly apply the preconditioner to random vectors and check
+  // that the preconditioned value is the same as one obtained by
+  // solving the linear system directly.
+  for (int i = 0; i < 5; ++i) {
+    CompressedRowSparseMatrix* lhs = inner_product_computer_->mutable_result();
+    Vector rhs = Vector::Random(lhs->num_rows());
+    Vector expected(lhs->num_rows());
+    EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));
+
+    Vector actual(lhs->num_rows());
+    preconditioner_->RightMultiply(rhs.data(), actual.data());
+
+    Matrix eigen_lhs;
+    lhs->ToDenseMatrix(&eigen_lhs);
+    EXPECT_NEAR((actual - expected).norm() / actual.norm(),
+                0.0,
+                std::numeric_limits<double>::epsilon() * 10)
+        << "\n"
+        << eigen_lhs << "\n"
+        << expected.transpose() << "\n"
+        << actual.transpose();
+  }
+}
+
+#ifndef CERES_NO_SUITESPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithSuiteSparse,
+                        SubsetPreconditionerTest,
+                        ::testing::Combine(::testing::Values(SUITE_SPARSE),
+                                           ::testing::Values(true, false)),
+                        ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_CXSPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithCXSparse,
+                        SubsetPreconditionerTest,
+                        ::testing::Combine(::testing::Values(CX_SPARSE),
+                                           ::testing::Values(true, false)),
+                        ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_EIGEN_SPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithEigenSparse,
+                        SubsetPreconditionerTest,
+                        ::testing::Combine(::testing::Values(EIGEN_SPARSE),
+                                           ::testing::Values(true, false)),
+                        ParamInfoToString);
+#endif
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/jni/Android.mk b/jni/Android.mk
index 8dbb593..4ddcdf8 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -195,6 +195,7 @@
                    $(CERES_SRC_PATH)/sparse_normal_cholesky_solver.cc \
                    $(CERES_SRC_PATH)/split.cc \
                    $(CERES_SRC_PATH)/stringprintf.cc \
+                   $(CERES_SRC_PATH)/subset_preconditioner.cc \
                    $(CERES_SRC_PATH)/suitesparse.cc \
                    $(CERES_SRC_PATH)/triplet_sparse_matrix.cc \
                    $(CERES_SRC_PATH)/trust_region_minimizer.cc \