diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 8a56ca0..052e7a0 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -70,6 +70,7 @@
     gradient_problem.cc
     gradient_problem_solver.cc
     implicit_schur_complement.cc
+    inner_product_computer.cc
     is_close.cc
     iterative_schur_complement_solver.cc
     levenberg_marquardt_strategy.cc
@@ -316,6 +317,7 @@
   ceres_test(graph_algorithms)
   ceres_test(householder_vector)
   ceres_test(implicit_schur_complement)
+  ceres_test(inner_product_computer)
   ceres_test(invert_psd_matrix)
   ceres_test(is_close)
   ceres_test(iterative_schur_complement_solver)
diff --git a/internal/ceres/inner_product_computer.cc b/internal/ceres/inner_product_computer.cc
new file mode 100644
index 0000000..052bd5a
--- /dev/null
+++ b/internal/ceres/inner_product_computer.cc
@@ -0,0 +1,373 @@
+// 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/inner_product_computer.h"
+
+#include <algorithm>
+
+namespace ceres {
+namespace internal {
+namespace {
+
+// Compute the product (in MATLAB notation)
+//
+// c(0:a_cols, 0:b_cols) = a' * b
+//
+// Where:
+//
+//  a is ab_rows x a_cols
+//  b is ab_rows x b_cols
+//  c is a_cos x c_col_stride
+//
+// a, b and c are row-major matrices.
+//
+// Performance note:
+// ----------------
+//
+// Technically this function is a repeat of a similarly named function
+// in small_blas.h but its performance is considerably better than
+// that of the version there due to the way it accesses memory.
+//
+// TODO(sameeragarwal): Measure and tune the performance of
+// small_blas.h based on the insights gained here.
+EIGEN_STRONG_INLINE void MatrixTransposeMatrixMultiply(const int ab_rows,
+                                                       const double* a,
+                                                       const int a_cols,
+                                                       const double* b,
+                                                       const int b_cols,
+                                                       double* c,
+                                                       int c_cols) {
+  // Compute c as the sum of ab_rows, rank 1 outer products of the
+  // corresponding rows of a and b.
+  for (int r = 0; r < ab_rows; ++r) {
+    double* c_r = c;
+    for (int i1 = 0; i1 < a_cols; ++i1) {
+      const double a_v = a[i1];
+      for (int i2 = 0; i2 < b_cols; ++i2) {
+        c_r[i2] += a_v * b[i2];
+      }
+      c_r += c_cols;
+    }
+    a += a_cols;
+    b += b_cols;
+  }
+}
+
+}  // namespace
+
+// Create the CompressedRowSparseMatrix matrix that will contain the
+// inner product.
+//
+// storage_type controls whether the result matrix contains the upper
+// or the lower triangular part of the product.
+//
+// num_nonzeros is the number of non-zeros in the result matrix.
+CompressedRowSparseMatrix* InnerProductComputer::CreateResultMatrix(
+    const CompressedRowSparseMatrix::StorageType storage_type,
+    const int num_nonzeros) {
+  CompressedRowSparseMatrix* matrix =
+      new CompressedRowSparseMatrix(m_.num_cols(), m_.num_cols(), num_nonzeros);
+  matrix->set_storage_type(storage_type);
+
+  const CompressedRowBlockStructure* bs = m_.block_structure();
+  const std::vector<Block>& blocks = bs->cols;
+  matrix->mutable_row_blocks()->resize(blocks.size());
+  matrix->mutable_col_blocks()->resize(blocks.size());
+  for (int i = 0; i < blocks.size(); ++i) {
+    (*(matrix->mutable_row_blocks()))[i] = blocks[i].size;
+    (*(matrix->mutable_col_blocks()))[i] = blocks[i].size;
+  }
+
+  return matrix;
+}
+
+// Given the set of product terms in the inner product, return the
+// total number of non-zeros in the result and for each row block of
+// the result matrix, compute the number of non-zeros in any one row
+// of the row block.
+int InnerProductComputer::ComputeNonzeros(
+    const std::vector<InnerProductComputer::ProductTerm>& product_terms,
+    std::vector<int>* row_nnz) {
+  const CompressedRowBlockStructure* bs = m_.block_structure();
+  const std::vector<Block>& blocks = bs->cols;
+
+  row_nnz->resize(blocks.size());
+  std::fill(row_nnz->begin(), row_nnz->end(), 0);
+
+  // First product term.
+  (*row_nnz)[product_terms[0].row] = blocks[product_terms[0].col].size;
+  int num_nonzeros =
+      blocks[product_terms[0].row].size * blocks[product_terms[0].col].size;
+
+  // Remaining product terms.
+  for (int i = 1; i < product_terms.size(); ++i) {
+    const ProductTerm& previous = product_terms[i - 1];
+    const ProductTerm& current = product_terms[i];
+
+    // Each (row, col) block counts only once.
+    // This check depends on product sorted on (row, col).
+    if (current.row != previous.row || current.col != previous.col) {
+      (*row_nnz)[current.row] += blocks[current.col].size;
+      num_nonzeros += blocks[current.row].size * blocks[current.col].size;
+    }
+  }
+
+  return num_nonzeros;
+}
+
+InnerProductComputer::InnerProductComputer(const BlockSparseMatrix& m,
+                                           const int start_row_block,
+                                           const int end_row_block)
+    : m_(m), start_row_block_(start_row_block), end_row_block_(end_row_block) {}
+
+// Compute the sparsity structure of the product m.transpose() * m
+// and create a CompressedRowSparseMatrix corresponding to it.
+//
+// Also compute the "program" vector, which for every term in the
+// block outer product provides the information for the entry in the
+// values array of the result matrix where it should be accumulated.
+//
+// Since the entries of the program are the same for rows with the
+// same sparsity structure, the program only stores the result for one
+// row per row block. The Compute function reuses this information for
+// each row in the row block.
+//
+// product_storage_type controls the form of the output matrix. It
+// can be LOWER_TRIANGULAR or UPPER_TRIANGULAR.
+InnerProductComputer* InnerProductComputer::Create(
+    const BlockSparseMatrix& m,
+    CompressedRowSparseMatrix::StorageType product_storage_type) {
+  return InnerProductComputer::Create(
+      m, 0, m.block_structure()->rows.size(), product_storage_type);
+}
+
+InnerProductComputer* InnerProductComputer::Create(
+    const BlockSparseMatrix& m,
+    const int start_row_block,
+    const int end_row_block,
+    CompressedRowSparseMatrix::StorageType product_storage_type) {
+  CHECK(product_storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR ||
+        product_storage_type == CompressedRowSparseMatrix::UPPER_TRIANGULAR);
+  CHECK_GT(m.num_nonzeros(), 0)
+      << "Congratulations, you found a bug in Ceres. Please report it.";
+  InnerProductComputer* inner_product_computer =
+      new InnerProductComputer(m, start_row_block, end_row_block);
+  inner_product_computer->Init(product_storage_type);
+  return inner_product_computer;
+}
+
+void InnerProductComputer::Init(
+    const CompressedRowSparseMatrix::StorageType product_storage_type) {
+  std::vector<InnerProductComputer::ProductTerm> product_terms;
+  const CompressedRowBlockStructure* bs = m_.block_structure();
+
+  // Give input matrix m in Block Sparse format
+  //     (row_block, col_block)
+  // represent each block multiplication
+  //     (row_block, col_block1)' X (row_block, col_block2)
+  // by its product term:
+  //     (col_block1, col_block2, index)
+  for (int row_block = start_row_block_; row_block < end_row_block_;
+       ++row_block) {
+    const CompressedRow& row = bs->rows[row_block];
+    for (int c1 = 0; c1 < row.cells.size(); ++c1) {
+      const Cell& cell1 = row.cells[c1];
+      int c2_begin, c2_end;
+      if (product_storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+        c2_begin = 0;
+        c2_end = c1 + 1;
+      } else {
+        c2_begin = c1;
+        c2_end = row.cells.size();
+      }
+
+      for (int c2 = c2_begin; c2 < c2_end; ++c2) {
+        const Cell& cell2 = row.cells[c2];
+        product_terms.push_back(InnerProductComputer::ProductTerm(
+            cell1.block_id, cell2.block_id, product_terms.size()));
+      }
+    }
+  }
+
+  std::sort(product_terms.begin(), product_terms.end());
+  ComputeOffsetsAndCreateResultMatrix(product_storage_type, product_terms);
+}
+
+void InnerProductComputer::ComputeOffsetsAndCreateResultMatrix(
+    const CompressedRowSparseMatrix::StorageType product_storage_type,
+    const std::vector<InnerProductComputer::ProductTerm>& product_terms) {
+  const int num_cols = m_.num_cols();
+  const std::vector<Block>& col_blocks = m_.block_structure()->cols;
+
+  std::vector<int> row_block_nnz;
+  const int num_nonzeros = ComputeNonzeros(product_terms, &row_block_nnz);
+
+  result_.reset(CreateResultMatrix(product_storage_type, num_nonzeros));
+
+  // Populate the row non-zero counts in the result matrix.
+  int* crsm_rows = result_->mutable_rows();
+  crsm_rows[0] = 0;
+  for (int i = 0; i < col_blocks.size(); ++i) {
+    for (int j = 0; j < col_blocks[i].size; ++j, ++crsm_rows) {
+      *(crsm_rows + 1) = *crsm_rows + row_block_nnz[i];
+    }
+  }
+
+  // The following macro FILL_CRSM_COL_BLOCK is key to understanding
+  // how this class works.
+  //
+  // It does two things.
+  //
+  // Sets the value for the current term in the result_offsets_ array
+  // and populates the cols array of the result matrix.
+  //
+  // row_block and col_block as the names imply, refer to the row and
+  // column blocks of the current term.
+  //
+  // row_nnz is the number of nonzeros in the result_matrix at the
+  // beginning of the first row of row_block.
+  //
+  // col_nnz is the number of nonzeros in the first row of the row
+  // block that occur before the current column block, i.e. this is
+  // sum of the sizes of all the column blocks in this row block that
+  // came before this column block.
+  //
+  // Given these two numbers and the total number of nonzeros in this
+  // row (nnz_in_row), we can now populate the cols array as follows:
+  //
+  // nnz + j * nnz_in_row is the beginning of the j^th row.
+  //
+  // nnz + j * nnz_in_row + col_nnz is the beginning of the column
+  // block in the j^th row.
+  //
+  // nnz + j * nnz_in_row + col_nnz + k is then the j^th row and the
+  // k^th column of the product block, whose value is
+  //
+  // col_blocks[col_block].position + k, which is the column number of
+  // the k^th column of the current column block.
+#define FILL_CRSM_COL_BLOCK                                \
+  const int row_block = current->row;                      \
+  const int col_block = current->col;                      \
+  const int nnz_in_row = row_block_nnz[row_block];         \
+  int* crsm_cols = result_->mutable_cols();                \
+  result_offsets_[current->index] = nnz + col_nnz;         \
+  for (int j = 0; j < col_blocks[row_block].size; ++j) {   \
+    for (int k = 0; k < col_blocks[col_block].size; ++k) { \
+      crsm_cols[nnz + j * nnz_in_row + col_nnz + k] =      \
+          col_blocks[col_block].position + k;              \
+    }                                                      \
+  }
+
+  result_offsets_.resize(product_terms.size());
+  int col_nnz = 0;
+  int nnz = 0;
+
+  // Process the first term.
+  const InnerProductComputer::ProductTerm* current = &product_terms[0];
+  FILL_CRSM_COL_BLOCK;
+
+  // Process the rest of the terms.
+  for (int i = 1; i < product_terms.size(); ++i) {
+    current = &product_terms[i];
+    const InnerProductComputer::ProductTerm* previous = &product_terms[i - 1];
+
+    // If the current term is the same as the previous term, then it
+    // stores its product at the same location as the previous term.
+    if (previous->row == current->row && previous->col == current->col) {
+      result_offsets_[current->index] = result_offsets_[previous->index];
+      continue;
+    }
+
+    if (previous->row == current->row) {
+      // if the current and previous terms are in the same row block,
+      // then they differ in the column block, in which case advance
+      // col_nnz by the column size of the prevous term.
+      col_nnz += col_blocks[previous->col].size;
+    } else {
+      // If we have moved to a new row-block , then col_nnz is zero,
+      // and nnz is set to the beginning of the row block.
+      col_nnz = 0;
+      nnz += row_block_nnz[previous->row] * col_blocks[previous->row].size;
+    }
+
+    FILL_CRSM_COL_BLOCK;
+  }
+}
+
+// Use the results_offsets_ array to numerically compute the product
+// m' * m and store it in result_.
+//
+// TODO(sameeragarwal): Multithreading support.
+void InnerProductComputer::Compute() {
+  const double* m_values = m_.values();
+  const CompressedRowBlockStructure* bs = m_.block_structure();
+
+  const CompressedRowSparseMatrix::StorageType storage_type =
+      result_->storage_type();
+  result_->SetZero();
+  double* values = result_->mutable_values();
+  const int* rows = result_->rows();
+  int cursor = 0;
+
+  // Iterate row blocks.
+  for (int r = start_row_block_; r < end_row_block_; ++r) {
+    const CompressedRow& m_row = bs->rows[r];
+    for (int c1 = 0; c1 < m_row.cells.size(); ++c1) {
+      const Cell& cell1 = m_row.cells[c1];
+      const int c1_size = bs->cols[cell1.block_id].size;
+      const int row_nnz = rows[bs->cols[cell1.block_id].position + 1] -
+          rows[bs->cols[cell1.block_id].position];
+
+      int c2_begin, c2_end;
+      if (storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+        c2_begin = 0;
+        c2_end = c1 + 1;
+      } else {
+        c2_begin = c1;
+        c2_end = m_row.cells.size();
+      }
+
+      for (int c2 = c2_begin; c2 < c2_end; ++c2, ++cursor) {
+        const Cell& cell2 = m_row.cells[c2];
+        const int c2_size = bs->cols[cell2.block_id].size;
+        MatrixTransposeMatrixMultiply(m_row.block.size,
+                                      m_values + cell1.position, c1_size,
+                                      m_values + cell2.position, c2_size,
+                                      values + result_offsets_[cursor],
+                                      row_nnz);
+      }
+    }
+  }
+
+  CHECK_EQ(cursor, result_offsets_.size());
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/inner_product_computer.h b/internal/ceres/inner_product_computer.h
new file mode 100644
index 0000000..d206707
--- /dev/null
+++ b/internal/ceres/inner_product_computer.h
@@ -0,0 +1,157 @@
+// 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_INNER_PRODUCT_COMPUTER_H_
+#define CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_
+
+#include <vector>
+
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/internal/scoped_ptr.h"
+
+namespace ceres {
+namespace internal {
+
+// This class is used to repeatedly compute the inner product
+//
+//   result = m' * m
+//
+// where the sparsity structure of m remains constant across calls.
+//
+// Upon creation, the class computes and caches information needed to
+// compute v, and then uses it to efficiently compute the product
+// every time InnerProductComputer::Compute is called.
+//
+// See sparse_normal_cholesky_solver.cc for example usage.
+//
+// Note that the result matrix is a block upper or lower-triangular
+// matrix, i.e., it will contain entries in the upper or lower
+// triangular part of the matrix corresponding to the block that occur
+// along its diagonal.
+//
+// This is not a problem as sparse linear algebra libraries can ignore
+// these entries with ease and the space used is minimal/linear in the
+// size of the matrices.
+class InnerProductComputer {
+ public:
+  // Factory
+  //
+  // m is the input matrix
+  //
+  // Since m' * m is a symmetric matrix, we only compute half of the
+  // matrix and the value of storage_type which must be
+  // UPPER_TRIANGULAR or LOWER_TRIANGULAR determines which half is
+  // computed.
+  //
+  // The user must ensure that the matrix m is valid for the life time
+  // of this object.
+  static InnerProductComputer* Create(
+      const BlockSparseMatrix& m,
+      CompressedRowSparseMatrix::StorageType storage_type);
+
+  // This factory method allows the user control over range of row
+  // blocks of m that should be used to compute the inner product.
+  //
+  // a = m(start_row_block : end_row_block, :);
+  // result = a' * a;
+  static InnerProductComputer* Create(
+      const BlockSparseMatrix& m,
+      int start_row_block,
+      int end_row_block,
+      CompressedRowSparseMatrix::StorageType storage_type);
+
+  // Update result_ to be numerically equal to m' * m.
+  void Compute();
+
+  // Accessors for the result containing the inner product.
+  //
+  // Compute must be called before accessing this result for
+  // the first time.
+  const CompressedRowSparseMatrix& result() const { return *result_; }
+  CompressedRowSparseMatrix* mutable_result() const { return result_.get(); }
+
+ private:
+  // A ProductTerm is a term in the block inner product of a matrix
+  // with itself.
+  struct ProductTerm {
+    ProductTerm(const int row, const int col, const int index)
+        : row(row), col(col), index(index) {}
+
+    bool operator<(const ProductTerm& right) const {
+      if (row == right.row) {
+        if (col == right.col) {
+          return index < right.index;
+        }
+        return col < right.col;
+      }
+      return row < right.row;
+    }
+
+    int row;
+    int col;
+    int index;
+  };
+
+  InnerProductComputer(const BlockSparseMatrix& m,
+                       int start_row_block,
+                       int end_row_block);
+
+  void Init(CompressedRowSparseMatrix::StorageType storage_type);
+
+  CompressedRowSparseMatrix* CreateResultMatrix(
+      const CompressedRowSparseMatrix::StorageType storage_type,
+      int num_nonzeros);
+
+  int ComputeNonzeros(const std::vector<ProductTerm>& product_terms,
+                      std::vector<int>* row_block_nnz);
+
+  void ComputeOffsetsAndCreateResultMatrix(
+      const CompressedRowSparseMatrix::StorageType storage_type,
+      const std::vector<ProductTerm>& product_terms);
+
+  const BlockSparseMatrix& m_;
+  const int start_row_block_;
+  const int end_row_block_;
+  scoped_ptr<CompressedRowSparseMatrix> result_;
+
+  // For each term in the inner product, result_offsets_ contains the
+  // location in the values array of the result_ matrix where it
+  // should be stored.
+  //
+  // This is the principal look up table that allows this class to
+  // compute the inner product fast.
+  std::vector<int> result_offsets_;
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_
diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc
new file mode 100644
index 0000000..a4c81bc
--- /dev/null
+++ b/internal/ceres/inner_product_computer_test.cc
@@ -0,0 +1,236 @@
+// 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/inner_product_computer.h"
+
+#include <numeric>
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/random.h"
+#include "ceres/triplet_sparse_matrix.h"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+#include "Eigen/SparseCore"
+
+namespace ceres {
+namespace internal {
+
+template <enum Eigen::UpLoType T>
+void CompareTriangularPartOfMatrices(const Matrix& expected,
+                                     const Matrix& actual) {
+  EXPECT_EQ(actual.rows(), actual.cols());
+  EXPECT_EQ(expected.rows(), expected.cols());
+  EXPECT_EQ(actual.rows(), expected.rows());
+
+  const Matrix expected_t = expected.triangularView<T>();
+  const Matrix actual_t = actual.triangularView<T>();
+
+  // TODO(sameeragarwal): Foo
+  CHECK_LE((expected_t - actual_t).norm() / actual_t.norm(),
+           100 * std::numeric_limits<double>::epsilon())
+      << "expected: \n"
+      << expected_t << "\nactual: \n"
+      << actual_t;
+}
+
+#define COMPUTE_AND_COMPARE                                  \
+  {                                                          \
+    inner_product_computer->Compute();                       \
+    CompressedRowSparseMatrix* actual_product_crsm =         \
+        inner_product_computer->mutable_result();            \
+    Matrix actual_inner_product =                            \
+        Eigen::MappedSparseMatrix<double, Eigen::ColMajor>(  \
+            actual_product_crsm->num_rows(),                 \
+            actual_product_crsm->num_rows(),                 \
+            actual_product_crsm->num_nonzeros(),             \
+            actual_product_crsm->mutable_rows(),             \
+            actual_product_crsm->mutable_cols(),             \
+            actual_product_crsm->mutable_values());          \
+    (actual_product_crsm->storage_type() ==                  \
+     CompressedRowSparseMatrix::LOWER_TRIANGULAR)            \
+        ? CompareTriangularPartOfMatrices<Eigen::Upper>(     \
+              expected_inner_product, actual_inner_product)  \
+        : CompareTriangularPartOfMatrices<Eigen::Lower>(     \
+              expected_inner_product, actual_inner_product); \
+  }
+
+
+TEST(InnerProductComputer, NormalOperation) {
+  // "Randomly generated seed."
+  SetRandomState(29823);
+  const int kMaxNumRowBlocks = 10;
+  const int kMaxNumColBlocks = 10;
+  const int kNumTrials = 10;
+
+  // Create a random matrix, compute its outer product using Eigen and
+  // ComputeOuterProduct. Convert both matrices to dense matrices and
+  // compare their upper triangular parts.
+  for (int num_row_blocks = 1; num_row_blocks < kMaxNumRowBlocks;
+       ++num_row_blocks) {
+    for (int num_col_blocks = 1; num_col_blocks < kMaxNumColBlocks;
+         ++num_col_blocks) {
+      for (int trial = 0; trial < kNumTrials; ++trial) {
+        BlockSparseMatrix::RandomMatrixOptions options;
+        options.num_row_blocks = num_row_blocks;
+        options.num_col_blocks = num_col_blocks;
+        options.min_row_block_size = 1;
+        options.max_row_block_size = 5;
+        options.min_col_block_size = 1;
+        options.max_col_block_size = 10;
+        options.block_density = std::max(0.1, RandDouble());
+
+        VLOG(2) << "num row blocks: " << options.num_row_blocks;
+        VLOG(2) << "num col blocks: " << options.num_col_blocks;
+        VLOG(2) << "min row block size: " << options.min_row_block_size;
+        VLOG(2) << "max row block size: " << options.max_row_block_size;
+        VLOG(2) << "min col block size: " << options.min_col_block_size;
+        VLOG(2) << "max col block size: " << options.max_col_block_size;
+        VLOG(2) << "block density: " << options.block_density;
+
+        scoped_ptr<BlockSparseMatrix> random_matrix(
+            BlockSparseMatrix::CreateRandomMatrix(options));
+
+        TripletSparseMatrix tsm(random_matrix->num_rows(),
+                                random_matrix->num_cols(),
+                                random_matrix->num_nonzeros());
+        random_matrix->ToTripletSparseMatrix(&tsm);
+        std::vector<Eigen::Triplet<double> > triplets;
+        for (int i = 0; i < tsm.num_nonzeros(); ++i) {
+          triplets.push_back(Eigen::Triplet<double>(
+              tsm.rows()[i], tsm.cols()[i], tsm.values()[i]));
+        }
+        Eigen::SparseMatrix<double> eigen_random_matrix(
+            random_matrix->num_rows(), random_matrix->num_cols());
+        eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end());
+        Matrix expected_inner_product =
+            eigen_random_matrix.transpose() * eigen_random_matrix;
+
+        scoped_ptr<InnerProductComputer> inner_product_computer;
+
+        inner_product_computer.reset(InnerProductComputer::Create(
+            *random_matrix, CompressedRowSparseMatrix::LOWER_TRIANGULAR));
+        COMPUTE_AND_COMPARE;
+        inner_product_computer.reset(InnerProductComputer::Create(
+            *random_matrix, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+        COMPUTE_AND_COMPARE;
+
+      }
+    }
+  }
+}
+
+
+TEST(InnerProductComputer, SubMatrix) {
+  // "Randomly generated seed."
+  SetRandomState(29823);
+  const int kNumRowBlocks = 10;
+  const int kNumColBlocks = 20;
+  const int kNumTrials = 5;
+
+  // Create a random matrix, compute its outer product using Eigen and
+  // ComputeInnerProductComputer. Convert both matrices to dense matrices and
+  // compare their upper triangular parts.
+  for (int trial = 0; trial < kNumTrials; ++trial) {
+    BlockSparseMatrix::RandomMatrixOptions options;
+    options.num_row_blocks = kNumRowBlocks;
+    options.num_col_blocks = kNumColBlocks;
+    options.min_row_block_size = 1;
+    options.max_row_block_size = 5;
+    options.min_col_block_size = 1;
+    options.max_col_block_size = 10;
+    options.block_density = std::max(0.1, RandDouble());
+
+    VLOG(2) << "num row blocks: " << options.num_row_blocks;
+    VLOG(2) << "num col blocks: " << options.num_col_blocks;
+    VLOG(2) << "min row block size: " << options.min_row_block_size;
+    VLOG(2) << "max row block size: " << options.max_row_block_size;
+    VLOG(2) << "min col block size: " << options.min_col_block_size;
+    VLOG(2) << "max col block size: " << options.max_col_block_size;
+    VLOG(2) << "block density: " << options.block_density;
+
+    scoped_ptr<BlockSparseMatrix> random_matrix(
+        BlockSparseMatrix::CreateRandomMatrix(options));
+
+    const std::vector<CompressedRow>& row_blocks =
+        random_matrix->block_structure()->rows;
+    const int num_row_blocks = row_blocks.size();
+
+    for (int start_row_block = 0; start_row_block < num_row_blocks - 1;
+         ++start_row_block) {
+      for (int end_row_block = start_row_block + 1;
+           end_row_block < num_row_blocks;
+           ++end_row_block) {
+        const int start_row = row_blocks[start_row_block].block.position;
+        const int end_row = row_blocks[end_row_block].block.position;
+
+        TripletSparseMatrix tsm(random_matrix->num_rows(),
+                                random_matrix->num_cols(),
+                                random_matrix->num_nonzeros());
+        random_matrix->ToTripletSparseMatrix(&tsm);
+        std::vector<Eigen::Triplet<double> > triplets;
+        for (int i = 0; i < tsm.num_nonzeros(); ++i) {
+          if (tsm.rows()[i] >= start_row && tsm.rows()[i] < end_row) {
+            triplets.push_back(Eigen::Triplet<double>(
+                tsm.rows()[i], tsm.cols()[i], tsm.values()[i]));
+          }
+        }
+
+        Eigen::SparseMatrix<double> eigen_random_matrix(
+            random_matrix->num_rows(), random_matrix->num_cols());
+        eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end());
+
+        Matrix expected_inner_product =
+            eigen_random_matrix.transpose() * eigen_random_matrix;
+
+        scoped_ptr<InnerProductComputer> inner_product_computer;
+        inner_product_computer.reset(InnerProductComputer::Create(
+            *random_matrix,
+            start_row_block,
+            end_row_block,
+            CompressedRowSparseMatrix::LOWER_TRIANGULAR));
+        COMPUTE_AND_COMPARE;
+        inner_product_computer.reset(InnerProductComputer::Create(
+            *random_matrix,
+            start_row_block,
+            end_row_block,
+            CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+        COMPUTE_AND_COMPARE;
+
+      }
+    }
+  }
+}
+
+#undef COMPUTE_AND_COMPARE
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/jni/Android.mk b/jni/Android.mk
index 20b484d..551cd45 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -154,6 +154,7 @@
                    $(CERES_SRC_PATH)/gradient_problem_solver.cc \
                    $(CERES_SRC_PATH)/is_close.cc \
                    $(CERES_SRC_PATH)/implicit_schur_complement.cc \
+                   $(CERES_SRC_PATH)/inner_product_computer.cc \
                    $(CERES_SRC_PATH)/iterative_schur_complement_solver.cc \
                    $(CERES_SRC_PATH)/lapack.cc \
                    $(CERES_SRC_PATH)/levenberg_marquardt_strategy.cc \