diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 439270e..34f31ab 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -125,7 +125,7 @@
 
   // Find the cumulative sum of the row counts.
   for (int i = 1; i < num_rows_ + 1; ++i) {
-    rows_[i] += rows_[i-1];
+    rows_[i] += rows_[i - 1];
   }
 
   CHECK_EQ(num_nonzeros(), m.num_nonzeros());
@@ -217,8 +217,8 @@
   num_rows_ -= delta_rows;
   rows_.resize(num_rows_ + 1);
 
-  // Walk the list of row blocks untill we reach the new number of
-  // rows and then drop the rest of the row blocks.
+  // Walk the list of row blocks until we reach the new number of rows
+  // and the drop the rest of the row blocks.
   int num_row_blocks = 0;
   int num_rows = 0;
   while (num_row_blocks < row_blocks_.size() && num_rows < num_rows_) {
@@ -380,6 +380,155 @@
   return transpose;
 }
 
+namespace {
+// A ProductTerm is a term in the outer 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;
+};
+
+CompressedRowSparseMatrix*
+CompressAndFillProgram(const int num_rows,
+                       const int num_cols,
+                       const vector<ProductTerm>& product,
+                       vector<int>* program) {
+  CHECK_GT(product.size(), 0);
+
+  // Count the number of unique product term, which in turn is the
+  // number of non-zeros in the outer product.
+  int num_nonzeros = 1;
+  for (int i = 1; i < product.size(); ++i) {
+    if (product[i].row != product[i - 1].row ||
+        product[i].col != product[i - 1].col) {
+      ++num_nonzeros;
+    }
+  }
+
+  CompressedRowSparseMatrix* matrix =
+      new CompressedRowSparseMatrix(num_rows, num_cols, num_nonzeros);
+
+  int* crsm_rows = matrix->mutable_rows();
+  std::fill(crsm_rows, crsm_rows + num_rows + 1, 0);
+  int* crsm_cols = matrix->mutable_cols();
+  std::fill(crsm_cols, crsm_cols + num_nonzeros, 0);
+
+  CHECK_NOTNULL(program)->clear();
+  program->resize(product.size());
+
+  // Iterate over the sorted product terms. This means each row is
+  // filled one at a time, and we are able to assign a position in the
+  // values array to each term.
+  //
+  // If terms repeat, i.e., they contribute to the same entry in the
+  // result matrix), then they do not affect the sparsity structure of
+  // the result matrix.
+  int nnz = 0;
+  crsm_cols[0] = product[0].col;
+  crsm_rows[product[0].row + 1]++;
+  (*program)[product[0].index] = nnz;
+  for (int i = 1; i < product.size(); ++i) {
+    const ProductTerm& previous = product[i - 1];
+    const ProductTerm& current = product[i];
+
+    // Sparsity structure is updated only if the term is not a repeat.
+    if (previous.row != current.row || previous.col != current.col) {
+      crsm_cols[++nnz] = current.col;
+      crsm_rows[current.row + 1]++;
+    }
+
+    // All terms get assigned the position in the values array where
+    // their value is accumulated.
+    (*program)[current.index] = nnz;
+  }
+
+  for (int i = 1; i < num_rows + 1; ++i) {
+    crsm_rows[i] += crsm_rows[i - 1];
+  }
+
+  return matrix;
+}
+
+}  // namespace
+
+CompressedRowSparseMatrix*
+CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
+      const CompressedRowSparseMatrix& m,
+      vector<int>* program) {
+  CHECK_NOTNULL(program)->clear();
+  CHECK_GT(m.num_nonzeros(), 0) << "Congratulations, "
+                                << "you found a bug in Ceres. Please report it.";
+
+  vector<ProductTerm> product;
+  const vector<int>& row_blocks = m.row_blocks();
+  int row_block_begin = 0;
+  // Iterate over row blocks
+  for (int row_block = 0; row_block < row_blocks.size(); ++row_block) {
+    const int row_block_end = row_block_begin + row_blocks[row_block];
+    // Compute the outer product terms for just one row per row block.
+    const int r = row_block_begin;
+    // Compute the lower triangular part of the product.
+    for (int idx1 = m.rows()[r]; idx1 < m.rows()[r + 1]; ++idx1) {
+      for (int idx2 = m.rows()[r]; idx2 <= idx1; ++idx2) {
+        product.push_back(ProductTerm(m.cols()[idx1], m.cols()[idx2], product.size()));
+      }
+    }
+    row_block_begin = row_block_end;
+  }
+  CHECK_EQ(row_block_begin, m.num_rows());
+  sort(product.begin(), product.end());
+  return CompressAndFillProgram(m.num_cols(), m.num_cols(), product, program);
+}
+
+void CompressedRowSparseMatrix::ComputeOuterProduct(
+    const CompressedRowSparseMatrix& m,
+    const vector<int>& program,
+    CompressedRowSparseMatrix* result) {
+  result->SetZero();
+  double* values = result->mutable_values();
+  const vector<int>& row_blocks = m.row_blocks();
+
+  int cursor = 0;
+  int row_block_begin = 0;
+  const double* m_values = m.values();
+  const int* m_rows = m.rows();
+  // Iterate over row blocks.
+  for (int row_block = 0; row_block < row_blocks.size(); ++row_block) {
+    const int row_block_end = row_block_begin + row_blocks[row_block];
+    const int saved_cursor = cursor;
+    for (int r = row_block_begin; r < row_block_end; ++r) {
+      // Reuse the program segment for each row in this row block.
+      cursor = saved_cursor;
+      const int row_begin = m_rows[r];
+      const int row_end = m_rows[r + 1];
+      for (int idx1 = row_begin; idx1 < row_end; ++idx1) {
+        const double v1 =  m_values[idx1];
+        for (int idx2 = row_begin; idx2 <= idx1; ++idx2, ++cursor) {
+          values[program[cursor]] += v1 * m_values[idx2];
+        }
+      }
+    }
+    row_block_begin = row_block_end;
+  }
+
+  CHECK_EQ(row_block_begin, m.num_rows());
+  CHECK_EQ(cursor, program.size());
+}
 
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index c5721eb..06b8689 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -128,6 +128,32 @@
       const double* diagonal,
       const vector<int>& blocks);
 
+  // Compute the sparsity structure of the product m.transpose() * m
+  // and create a CompressedRowSparseMatrix corresponding to it.
+  //
+  // Also compute a "program" vector, which for every term in the
+  // outer product points to the entry in the values array of the
+  // result matrix where it should be accumulated.
+  //
+  // This program is used by the ComputeOuterProduct function below to
+  // compute the outer product.
+  //
+  // 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 ComputeOuterProduct function reuses
+  // this information for each row in the row block.
+  static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram(
+      const CompressedRowSparseMatrix& m,
+      vector<int>* program);
+
+  // Compute the values array for the expression m.transpose() * m,
+  // where the matrix used to store the result and a program have been
+  // created using the CreateOuterProductMatrixAndProgram function
+  // above.
+  static void ComputeOuterProduct(const CompressedRowSparseMatrix& m,
+                                  const vector<int>& program,
+                                  CompressedRowSparseMatrix* result);
+
  private:
   int num_rows_;
   int num_cols_;
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index 4d20083..ea737b8 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -30,11 +30,14 @@
 
 #include "ceres/compressed_row_sparse_matrix.h"
 
+#include <numeric>
 #include "ceres/casts.h"
 #include "ceres/crs_matrix.h"
+#include "ceres/cxsparse.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/linear_least_squares_problems.h"
+#include "ceres/random.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "glog/logging.h"
 #include "gtest/gtest.h"
@@ -381,5 +384,170 @@
   EXPECT_NEAR((dense_matrix - dense_transpose.transpose()).norm(), 0.0, 1e-14);
 }
 
+#ifndef CERES_NO_CXSPARSE
+
+struct RandomMatrixOptions {
+  int num_row_blocks;
+  int min_row_block_size;
+  int max_row_block_size;
+  int num_col_blocks;
+  int min_col_block_size;
+  int max_col_block_size;
+  double block_density;
+};
+
+CompressedRowSparseMatrix* CreateRandomCompressedRowSparseMatrix(
+    const RandomMatrixOptions& options) {
+  vector<int> row_blocks;
+  for (int i = 0; i < options.num_row_blocks; ++i) {
+    const int delta_block_size =
+        Uniform(options.max_row_block_size - options.min_row_block_size);
+    row_blocks.push_back(options.min_row_block_size + delta_block_size);
+  }
+
+  vector<int> col_blocks;
+  for (int i = 0; i < options.num_col_blocks; ++i) {
+    const int delta_block_size =
+        Uniform(options.max_col_block_size - options.min_col_block_size);
+    col_blocks.push_back(options.min_col_block_size + delta_block_size);
+  }
+
+  vector<int> rows;
+  vector<int> cols;
+  vector<double> values;
+
+  while (values.size() == 0) {
+    int row_block_begin = 0;
+    for (int r = 0; r < options.num_row_blocks; ++r) {
+      int col_block_begin = 0;
+      for (int c = 0; c < options.num_col_blocks; ++c) {
+        if (RandDouble() <= options.block_density) {
+          for (int i = 0; i < row_blocks[r]; ++i) {
+            for (int j = 0; j < col_blocks[c]; ++j) {
+              rows.push_back(row_block_begin + i);
+              cols.push_back(col_block_begin + j);
+              values.push_back(RandNormal());
+            }
+          }
+        }
+        col_block_begin += col_blocks[c];
+      }
+      row_block_begin += row_blocks[r];
+    }
+  }
+
+  const int num_rows = std::accumulate(row_blocks.begin(), row_blocks.end(), 0);
+  const int num_cols = std::accumulate(col_blocks.begin(), col_blocks.end(), 0);
+  const int num_nonzeros = values.size();
+
+  TripletSparseMatrix tsm(num_rows, num_cols, num_nonzeros);
+  std::copy(rows.begin(), rows.end(), tsm.mutable_rows());
+  std::copy(cols.begin(), cols.end(), tsm.mutable_cols());
+  std::copy(values.begin(), values.end(), tsm.mutable_values());
+  tsm.set_num_nonzeros(num_nonzeros);
+  CompressedRowSparseMatrix* matrix = new CompressedRowSparseMatrix(tsm);
+  (*matrix->mutable_row_blocks())  = row_blocks;
+  (*matrix->mutable_col_blocks())  = col_blocks;
+  return matrix;
+}
+
+void ToDenseMatrix(const cs_di* matrix, Matrix* dense_matrix) {
+  dense_matrix->resize(matrix->m, matrix->n);
+  dense_matrix->setZero();
+
+  for (int c = 0; c < matrix->n; ++c) {
+   for (int idx = matrix->p[c]; idx < matrix->p[c + 1]; ++idx) {
+     const int r = matrix->i[idx];
+     (*dense_matrix)(r, c) = matrix->x[idx];
+   }
+ }
+}
+
+TEST(CompressedRowSparseMatrix, ComputeOuterProduct) {
+  // "Randomly generated seed."
+  SetRandomState(29823);
+  int kMaxNumRowBlocks = 10;
+  int kMaxNumColBlocks = 10;
+  int kNumTrials = 10;
+
+  CXSparse cxsparse;
+  const double kTolerance = 1e-18;
+
+  // Create a random matrix, compute its outer product using CXSParse
+  // and ComputeOuterProduct. Convert both matrices to dense matrices
+  // and compare their upper triangular parts. They should be within
+  // kTolerance of each other.
+  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) {
+
+
+        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<CompressedRowSparseMatrix> matrix(
+            CreateRandomCompressedRowSparseMatrix(options));
+
+        cs_di cs_matrix_transpose = cxsparse.CreateSparseMatrixTransposeView(matrix.get());
+        cs_di* cs_matrix = cxsparse.TransposeMatrix(&cs_matrix_transpose);
+        cs_di* expected_outer_product =
+            cxsparse.MatrixMatrixMultiply(&cs_matrix_transpose, cs_matrix);
+
+        vector<int> program;
+        scoped_ptr<CompressedRowSparseMatrix> outer_product(
+            CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
+                *matrix, &program));
+        CompressedRowSparseMatrix::ComputeOuterProduct(*matrix,
+                                                       program,
+                                                       outer_product.get());
+
+        cs_di actual_outer_product =
+            cxsparse.CreateSparseMatrixTransposeView(outer_product.get());
+
+        ASSERT_EQ(actual_outer_product.m, actual_outer_product.n);
+        ASSERT_EQ(expected_outer_product->m, expected_outer_product->n);
+        ASSERT_EQ(actual_outer_product.m, expected_outer_product->m);
+
+        Matrix actual_matrix;
+        Matrix expected_matrix;
+
+        ToDenseMatrix(expected_outer_product, &expected_matrix);
+        expected_matrix.triangularView<Eigen::StrictlyLower>().setZero();
+
+        ToDenseMatrix(&actual_outer_product, &actual_matrix);
+        const double diff_norm = (actual_matrix - expected_matrix).norm() / expected_matrix.norm();
+        ASSERT_NEAR(diff_norm, 0.0, kTolerance)
+            << "expected: \n"
+            << expected_matrix
+            << "\nactual: \n"
+            << actual_matrix;
+
+        cxsparse.Free(cs_matrix);
+        cxsparse.Free(expected_outer_product);
+      }
+    }
+  }
+}
+
+#endif  // CERES_NO_CXSPARSE
+
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index b2dba44..880adf7 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -77,37 +77,10 @@
     const double* b,
     const LinearSolver::PerSolveOptions& per_solve_options,
     double * x) {
-  switch (options_.sparse_linear_algebra_library_type) {
-    case SUITE_SPARSE:
-      return SolveImplUsingSuiteSparse(A, b, per_solve_options, x);
-    case CX_SPARSE:
-      return SolveImplUsingCXSparse(A, b, per_solve_options, x);
-    default:
-      LOG(FATAL) << "Unknown sparse linear algebra library : "
-                 << options_.sparse_linear_algebra_library_type;
-  }
-
-  LOG(FATAL) << "Unknown sparse linear algebra library : "
-             << options_.sparse_linear_algebra_library_type;
-  return LinearSolver::Summary();
-}
-
-#ifndef CERES_NO_CXSPARSE
-LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
-    CompressedRowSparseMatrix* A,
-    const double* b,
-    const LinearSolver::PerSolveOptions& per_solve_options,
-    double * x) {
-  EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve");
-
-  LinearSolver::Summary summary;
-  summary.num_iterations = 1;
-  summary.termination_type = LINEAR_SOLVER_SUCCESS;
-  summary.message = "Success.";
 
   const int num_cols = A->num_cols();
-  Vector Atb = Vector::Zero(num_cols);
-  A->LeftMultiply(b, Atb.data());
+  VectorRef(x, num_cols).setZero();
+  A->LeftMultiply(b, x);
 
   if (per_solve_options.D != NULL) {
     // Temporarily append a diagonal block to the A matrix, but undo
@@ -123,10 +96,37 @@
     A->AppendRows(*regularizer);
   }
 
-  VectorRef(x, num_cols).setZero();
+  LinearSolver::Summary summary;
+  switch (options_.sparse_linear_algebra_library_type) {
+    case SUITE_SPARSE:
+      summary = SolveImplUsingSuiteSparse(A, per_solve_options, x);
+      break;
+    case CX_SPARSE:
+      summary = SolveImplUsingCXSparse(A, per_solve_options, x);
+      break;
+    default:
+      LOG(FATAL) << "Unknown sparse linear algebra library : "
+                 << options_.sparse_linear_algebra_library_type;
+  }
 
-  // Wrap the augmented Jacobian in a compressed sparse column matrix.
-  cs_di At = cxsparse_.CreateSparseMatrixTransposeView(A);
+  if (per_solve_options.D != NULL) {
+    A->DeleteRows(num_cols);
+  }
+
+  return summary;
+}
+
+#ifndef CERES_NO_CXSPARSE
+LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
+    CompressedRowSparseMatrix* A,
+    const LinearSolver::PerSolveOptions& per_solve_options,
+    double * rhs_and_solution) {
+  EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve");
+
+  LinearSolver::Summary summary;
+  summary.num_iterations = 1;
+  summary.termination_type = LINEAR_SOLVER_SUCCESS;
+  summary.message = "Success.";
 
   // Compute the normal equations. J'J delta = J'f and solve them
   // using a sparse Cholesky factorization. Notice that when compared
@@ -135,13 +135,18 @@
   // factorized. CHOLMOD/SuiteSparse on the other hand can just work
   // off of Jt to compute the Cholesky factorization of the normal
   // equations.
-  cs_di* A2 = cxsparse_.TransposeMatrix(&At);
-  cs_di* AtA = cxsparse_.MatrixMatrixMultiply(&At, A2);
-
-  cxsparse_.Free(A2);
-  if (per_solve_options.D != NULL) {
-    A->DeleteRows(num_cols);
+  if (outer_product_.get() == NULL) {
+    outer_product_.reset(
+        CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
+            *A, &pattern_));
   }
+
+  CompressedRowSparseMatrix::ComputeOuterProduct(
+      *A, pattern_, outer_product_.get());
+  cs_di AtA_view =
+      cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
+  cs_di* AtA = &AtA_view;
+
   event_logger.AddEvent("Setup");
 
   // Compute symbolic factorization if not available.
@@ -160,23 +165,17 @@
     summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
     summary.message =
         "CXSparse failure. Unable to find symbolic factorization.";
-  } else if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) {
-    VectorRef(x, Atb.rows()) = Atb;
-  } else {
+  } else if (!cxsparse_.SolveCholesky(AtA, cxsparse_factor_, rhs_and_solution)) {
     summary.termination_type = LINEAR_SOLVER_FAILURE;
   }
   event_logger.AddEvent("Solve");
-
-  cxsparse_.Free(AtA);
-  event_logger.AddEvent("Teardown");
   return summary;
 }
 #else
 LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
     CompressedRowSparseMatrix* A,
-    const double* b,
     const LinearSolver::PerSolveOptions& per_solve_options,
-    double * x) {
+    double * rhs_and_solution) {
   LOG(FATAL) << "No CXSparse support in Ceres.";
 
   // Unreachable but MSVC does not know this.
@@ -187,9 +186,8 @@
 #ifndef CERES_NO_SUITESPARSE
 LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
     CompressedRowSparseMatrix* A,
-    const double* b,
     const LinearSolver::PerSolveOptions& per_solve_options,
-    double * x) {
+    double * rhs_and_solution) {
   EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve");
   LinearSolver::Summary summary;
   summary.termination_type = LINEAR_SOLVER_SUCCESS;
@@ -197,24 +195,6 @@
   summary.message = "Success.";
 
   const int num_cols = A->num_cols();
-  Vector Atb = Vector::Zero(num_cols);
-  A->LeftMultiply(b, Atb.data());
-
-  if (per_solve_options.D != NULL) {
-    // Temporarily append a diagonal block to the A matrix, but undo
-    // it before returning the matrix to the user.
-    scoped_ptr<CompressedRowSparseMatrix> regularizer;
-    if (A->col_blocks().size() > 0) {
-      regularizer.reset(CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
-                            per_solve_options.D, A->col_blocks()));
-    } else {
-      regularizer.reset(new CompressedRowSparseMatrix(
-                            per_solve_options.D, num_cols));
-    }
-    A->AppendRows(*regularizer);
-  }
-
-  VectorRef(x, num_cols).setZero();
   cholmod_sparse lhs = ss_.CreateSparseMatrixTransposeView(A);
   event_logger.AddEvent("Setup");
 
@@ -231,33 +211,23 @@
   event_logger.AddEvent("Analysis");
 
   if (factor_ == NULL) {
-    if (per_solve_options.D != NULL) {
-      A->DeleteRows(num_cols);
-    }
     summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
     return summary;
   }
 
   summary.termination_type = ss_.Cholesky(&lhs, factor_, &summary.message);
   if (summary.termination_type != LINEAR_SOLVER_SUCCESS) {
-    if (per_solve_options.D != NULL) {
-      A->DeleteRows(num_cols);
-    }
     return summary;
   }
 
-  cholmod_dense* rhs = ss_.CreateDenseVector(Atb.data(), num_cols, num_cols);
-  cholmod_dense* sol = ss_.Solve(factor_, rhs, &summary.message);
+  cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols);
+  cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message);
   event_logger.AddEvent("Solve");
 
   ss_.Free(rhs);
-  if (per_solve_options.D != NULL) {
-    A->DeleteRows(num_cols);
-  }
-
-  if (sol != NULL) {
-    memcpy(x, sol->x, num_cols * sizeof(*x));
-    ss_.Free(sol);
+  if (solution != NULL) {
+    memcpy(rhs_and_solution, solution->x, num_cols * sizeof(*rhs_and_solution));
+    ss_.Free(solution);
   } else {
     summary.termination_type = LINEAR_SOLVER_FAILURE;
   }
@@ -268,9 +238,8 @@
 #else
 LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
     CompressedRowSparseMatrix* A,
-    const double* b,
     const LinearSolver::PerSolveOptions& per_solve_options,
-    double * x) {
+    double * rhs_and_solution) {
   LOG(FATAL) << "No SuiteSparse support in Ceres.";
 
   // Unreachable but MSVC does not know this.
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index 61111b4..adca08a 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -62,16 +62,14 @@
 
   LinearSolver::Summary SolveImplUsingSuiteSparse(
       CompressedRowSparseMatrix* A,
-      const double* b,
       const LinearSolver::PerSolveOptions& options,
-      double* x);
+      double* rhs_and_solution);
 
   // Crashes if CSparse is not installed.
   LinearSolver::Summary SolveImplUsingCXSparse(
       CompressedRowSparseMatrix* A,
-      const double* b,
       const LinearSolver::PerSolveOptions& options,
-      double* x);
+      double* rhs_and_solution);
 
   SuiteSparse ss_;
   // Cached factorization
@@ -80,7 +78,8 @@
   CXSparse cxsparse_;
   // Cached factorization
   cs_dis* cxsparse_factor_;
-
+  scoped_ptr<CompressedRowSparseMatrix> outer_product_;
+  vector<int> pattern_;
   const LinearSolver::Options options_;
   CERES_DISALLOW_COPY_AND_ASSIGN(SparseNormalCholeskySolver);
 };
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc
index 949fdd1..624db68 100644
--- a/internal/ceres/unsymmetric_linear_solver_test.cc
+++ b/internal/ceres/unsymmetric_linear_solver_test.cc
@@ -57,7 +57,7 @@
   }
 
   void TestSolver(const LinearSolver::Options& options) {
-    scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
+
 
     LinearSolver::PerSolveOptions per_solve_options;
     LinearSolver::Summary unregularized_solve_summary;
@@ -84,13 +84,17 @@
     } else {
       LOG(FATAL) << "Unknown linear solver : " << options.type;
     }
+
     // Unregularized
+    scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
     unregularized_solve_summary =
         solver->Solve(transformed_A.get(),
                       b_.get(),
                       per_solve_options,
                       x_unregularized.data());
 
+    // Sparsity structure is changing, reset the solver.
+    solver.reset(LinearSolver::Create(options));
     // Regularized solution
     per_solve_options.D = D_.get();
     regularized_solve_summary =