diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 68af5f8..dceb3bc 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -109,6 +109,7 @@
                                                      int max_num_nonzeros) {
   num_rows_ = num_rows;
   num_cols_ = num_cols;
+  storage_type_ = UNSYMMETRIC;
   rows_.resize(num_rows + 1, 0);
   cols_.resize(max_num_nonzeros, 0);
   values_.resize(max_num_nonzeros, 0.0);
@@ -124,6 +125,7 @@
     const TripletSparseMatrix& m) {
   num_rows_ = m.num_rows();
   num_cols_ = m.num_cols();
+  storage_type_ = UNSYMMETRIC;
 
   rows_.resize(num_rows_ + 1, 0);
   cols_.resize(m.num_nonzeros(), 0);
@@ -168,6 +170,7 @@
 
   num_rows_ = num_rows;
   num_cols_ = num_rows;
+  storage_type_ = UNSYMMETRIC;
   rows_.resize(num_rows + 1);
   cols_.resize(num_rows);
   values_.resize(num_rows);
@@ -448,6 +451,20 @@
   CompressedRowSparseMatrix* transpose =
       new CompressedRowSparseMatrix(num_cols_, num_rows_, num_nonzeros());
 
+  switch (storage_type_) {
+    case UNSYMMETRIC:
+      transpose->set_storage_type(UNSYMMETRIC);
+      break;
+    case LOWER_TRIANGULAR:
+      transpose->set_storage_type(UPPER_TRIANGULAR);
+      break;
+    case UPPER_TRIANGULAR:
+      transpose->set_storage_type(LOWER_TRIANGULAR);
+      break;
+    default:
+      LOG(FATAL) << "Unknown storage type: " << storage_type_;
+  };
+
   TransposeForCompressedRowSparseStructure(num_rows(),
                                            num_cols(),
                                            num_nonzeros(),
@@ -524,6 +541,7 @@
 // the outerproduct matrix, which is used in outer product computation.
 CompressedRowSparseMatrix* CreateOuterProductMatrix(
     const int num_cols,
+    const CompressedRowSparseMatrix::StorageType storage_type,
     const vector<int>& blocks,
     const vector<ProductTerm>& product,
     vector<int>* row_nnz) {
@@ -546,6 +564,7 @@
 
   CompressedRowSparseMatrix* matrix =
       new CompressedRowSparseMatrix(num_cols, num_cols, num_nonzeros);
+  matrix->set_storage_type(storage_type);
 
   // Compute block offsets for outer product matrix, which is used
   // in ComputeOuterProduct.
@@ -561,6 +580,7 @@
 
 CompressedRowSparseMatrix* CompressAndFillProgram(
     const int num_cols,
+    const CompressedRowSparseMatrix::StorageType storage_type,
     const vector<int>& blocks,
     const vector<ProductTerm>& product,
     vector<int>* program) {
@@ -568,7 +588,7 @@
 
   vector<int> row_nnz;
   CompressedRowSparseMatrix* matrix =
-      CreateOuterProductMatrix(num_cols, blocks, product, &row_nnz);
+      CreateOuterProductMatrix(num_cols, storage_type, blocks, product, &row_nnz);
 
   const vector<int>& block_offsets = matrix->block_offsets();
 
@@ -679,7 +699,10 @@
 
 CompressedRowSparseMatrix*
 CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
-    const CompressedRowSparseMatrix& m, const int stype, vector<int>* program) {
+      const CompressedRowSparseMatrix& m,
+      const CompressedRowSparseMatrix::StorageType storage_type,
+      vector<int>* program) {
+  CHECK(storage_type ==  LOWER_TRIANGULAR || storage_type == UPPER_TRIANGULAR);
   CHECK_NOTNULL(program)->clear();
   CHECK_GT(m.num_nonzeros(), 0)
       << "Congratulations, you found a bug in Ceres. Please report it.";
@@ -701,7 +724,7 @@
   for (int row_block = 1; row_block < crsb_rows.size(); ++row_block) {
     for (int idx1 = crsb_rows[row_block - 1]; idx1 < crsb_rows[row_block];
          ++idx1) {
-      if (stype > 0) {  // Lower triangular matrix.
+      if (storage_type == LOWER_TRIANGULAR) {
         for (int idx2 = crsb_rows[row_block - 1]; idx2 <= idx1; ++idx2) {
           product.push_back(
               ProductTerm(crsb_cols[idx1], crsb_cols[idx2], product.size()));
@@ -716,7 +739,8 @@
   }
 
   sort(product.begin(), product.end());
-  return CompressAndFillProgram(m.num_cols(), col_blocks, product, program);
+  return CompressAndFillProgram(
+      m.num_cols(), storage_type, col_blocks, product, program);
 }
 
 // Give input matrix m in Compressed Row Sparse Block format
@@ -754,9 +778,10 @@
 // So there is no special handling for diagonal blocks.
 void CompressedRowSparseMatrix::ComputeOuterProduct(
     const CompressedRowSparseMatrix& m,
-    const int stype,
     const vector<int>& program,
     CompressedRowSparseMatrix* result) {
+  CHECK(result->storage_type() ==  LOWER_TRIANGULAR ||
+        result->storage_type() ==  UPPER_TRIANGULAR);
   result->SetZero();
   double* values = result->mutable_values();
   const int* rows = result->rows();
@@ -769,10 +794,11 @@
   const vector<int>& col_blocks = m.col_blocks();
   const vector<int>& crsb_rows = m.crsb_rows();
   const vector<int>& crsb_cols = m.crsb_cols();
-
+  const StorageType storage_type = result->storage_type();
 #define COL_BLOCK1 (crsb_cols[idx1])
 #define COL_BLOCK2 (crsb_cols[idx2])
 
+
   // Iterate row blocks.
   for (int row_block = 0, m_row_begin = 0; row_block < row_blocks.size();
        m_row_begin += row_blocks[row_block++]) {
@@ -790,8 +816,7 @@
       // outerproduct matrix compressed row.
       const int row_begin = block_offsets[COL_BLOCK1];
       const int row_nnz = rows[row_begin + 1] - rows[row_begin];
-
-      if (stype > 0) {  // Lower triangular matrix.
+      if (storage_type == LOWER_TRIANGULAR) {
         for (int idx2 = crsb_rows[row_block], m_col_nnz2 = 0; idx2 <= idx1;
              m_col_nnz2 += col_blocks[COL_BLOCK2], ++idx2, ++cursor) {
           int col_nnz = program[cursor];
@@ -805,7 +830,7 @@
                                      row_nnz,
                                      values + rows[row_begin] + col_nnz);
         }
-      } else {  // Upper triangular matrix.
+      } else {
         for (int idx2 = idx1, m_col_nnz2 = m_col_nnz1;
              idx2 < crsb_rows[row_block + 1];
              m_col_nnz2 += col_blocks[COL_BLOCK2], ++idx2, ++cursor) {
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index 0bc5d01..9759d13 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -48,6 +48,13 @@
 
 class CompressedRowSparseMatrix : public SparseMatrix {
  public:
+  enum StorageType {
+    UNSYMMETRIC,
+    LOWER_TRIANGULAR,
+    UPPER_TRIANGULAR
+  };
+
+
   // Build a matrix with the same content as the TripletSparseMatrix
   // m. TripletSparseMatrix objects are easier to construct
   // incrementally, so we use them to initialize SparseMatrix
@@ -102,6 +109,18 @@
 
   void ToCRSMatrix(CRSMatrix* matrix) const;
 
+  void SolveLowerTriangularInPlace(double* solution) const;
+  void SolveLowerTriangularTransposeInPlace(double* solution) const;
+
+  CompressedRowSparseMatrix* Transpose() const;
+
+  // Destructive array resizing method.
+  void SetMaxNumNonZeros(int num_nonzeros);
+
+  // Non-destructive array resizing method.
+  void set_num_rows(const int num_rows) { num_rows_ = num_rows; }
+  void set_num_cols(const int num_cols) { num_cols_ = num_cols; }
+
   // Low level access methods that expose the structure of the matrix.
   const int* cols() const { return &cols_[0]; }
   int* mutable_cols() { return &cols_[0]; }
@@ -109,6 +128,11 @@
   const int* rows() const { return &rows_[0]; }
   int* mutable_rows() { return &rows_[0]; }
 
+  const StorageType& storage_type() const { return storage_type_; }
+  void set_storage_type(const StorageType& storage_type) {
+    storage_type_ = storage_type;
+  }
+
   const std::vector<int>& row_blocks() const { return row_blocks_; }
   std::vector<int>* mutable_row_blocks() { return &row_blocks_; }
 
@@ -124,18 +148,6 @@
   const std::vector<int>& crsb_cols() const { return crsb_cols_; }
   std::vector<int>* mutable_crsb_cols() { return &crsb_cols_; }
 
-  // Destructive array resizing method.
-  void SetMaxNumNonZeros(int num_nonzeros);
-
-  // Non-destructive array resizing method.
-  void set_num_rows(const int num_rows) { num_rows_ = num_rows; }
-  void set_num_cols(const int num_cols) { num_cols_ = num_cols; }
-
-  void SolveLowerTriangularInPlace(double* solution) const;
-  void SolveLowerTriangularTransposeInPlace(double* solution) const;
-
-  CompressedRowSparseMatrix* Transpose() const;
-
   static CompressedRowSparseMatrix* CreateBlockDiagonalMatrix(
       const double* diagonal,
       const std::vector<int>& blocks);
@@ -155,12 +167,11 @@
   // one row per row block. The ComputeOuterProduct function reuses
   // this information for each row in the row block.
   //
-  // stype controls the result matrix in upper or lower triangular matrix
-  //    stype = 1:  lower triangular matrix
-  //    stype = -1: higher triangular matrix
+  // storage_type controls the form of the output matrix. It can be
+  // LOWER_TRIANGULAR or UPPER_TRIANGULAR.
   static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram(
       const CompressedRowSparseMatrix& m,
-      const int stype,
+      const StorageType storage_type,
       std::vector<int>* program);
 
   // Compute the values array for the expression m.transpose() * m,
@@ -168,7 +179,6 @@
   // created using the CreateOuterProductMatrixAndProgram function
   // above.
   static void ComputeOuterProduct(const CompressedRowSparseMatrix& m,
-                                  const int stype,
                                   const std::vector<int>& program,
                                   CompressedRowSparseMatrix* result);
 
@@ -178,6 +188,7 @@
   std::vector<int> rows_;
   std::vector<int> cols_;
   std::vector<double> values_;
+  StorageType storage_type_;
 
   // If the matrix has an underlying block structure, then it can also
   // carry with it row and column block sizes. This is auxilliary and
@@ -187,11 +198,11 @@
   std::vector<int> row_blocks_;
   std::vector<int> col_blocks_;
 
-  // For outerproduct matrix (J' * J), we pre-compute its block offsets
-  // information here for fast outerproduct computation in block unit.
-  // Since the outerproduct matrix is symmetric, we do not need to
-  // distinguish row or col block. In another word, this is the prefix
-  // sum of row_blocks_/col_blocks_.
+  // For outer product matrix (J' * J), we pre-compute its block
+  // offsets information here for fast outer product computation in
+  // block unit.  Since the outer product matrix is symmetric, we do
+  // not need to distinguish row or col block. In another word, this
+  // is the prefix sum of row_blocks_/col_blocks_.
   std::vector<int> block_offsets_;
 
   // If the matrix has an underlying block structure, then it can also
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index b6e732c..cc5ed08 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -494,16 +494,16 @@
         Matrix expected_outer_product =
             mapped_random_matrix.transpose() * mapped_random_matrix;
 
-        // Use compressed row lower triangular matrix.
-        const int stype = 1;
+        // Use compressed row lower triangular matrix, which will then
+        // get mapped to a compressed column upper triangular matrix.
         vector<int> program;
         scoped_ptr<CompressedRowSparseMatrix> outer_product(
             CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
-                *random_matrix, stype, &program));
-        CompressedRowSparseMatrix::ComputeOuterProduct(*random_matrix,
-                                                       stype,
-                                                       program,
-                                                       outer_product.get());
+                *random_matrix,
+                CompressedRowSparseMatrix::LOWER_TRIANGULAR,
+                &program));
+        CompressedRowSparseMatrix::ComputeOuterProduct(
+            *random_matrix, program, outer_product.get());
 
         Matrix actual_outer_product =
             Eigen::MappedSparseMatrix<double, Eigen::ColMajor>(
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 8401202..64fd080 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -275,21 +275,18 @@
                                &event_logger);
   }
 
-  // Compute outerproduct to compressed row lower triangular matrix.
-  // Eigen SimplicialLDLT default uses lower triangular part of matrix.
-  // This can change to upper triangular matrix if specifying
-  //    Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >
-  // with _UpLo = Upper.
-  const int stype = 1;
-
+  // Compute outer product as a compressed row lower triangular
+  // matrix, because after mapping to a column major matrix, this will
+  // become a compressed column upper triangular matrix. Which is the
+  // default that Eigen uses.
   if (outer_product_.get() == NULL) {
     outer_product_.reset(
         CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
-            *A, stype, &pattern_));
+            *A, CompressedRowSparseMatrix::LOWER_TRIANGULAR, &pattern_));
   }
 
   CompressedRowSparseMatrix::ComputeOuterProduct(
-      *A, stype, pattern_, outer_product_.get());
+      *A, pattern_, outer_product_.get());
 
   // Map to an upper triangular column major matrix.
   //
@@ -369,23 +366,21 @@
   summary.termination_type = LINEAR_SOLVER_SUCCESS;
   summary.message = "Success.";
 
-  // Compute outerproduct to compressed row lower triangular matrix.
-  // CXSparse Cholesky factorization uses lower triangular part of the matrix.
-  const int stype = 1;
 
-  // Compute the normal equations. J'J delta = J'f and solve them
-  // using a sparse Cholesky factorization. Notice that we explicitly
-  // compute the normal equations before they can be factorized.
+  // Compute outer product as a compressed row lower triangular
+  // matrix, which would mapped to a compressed column upper
+  // triangular matrix, which is the representation used by CXSparse's
+  // sparse Cholesky factorization.
   if (outer_product_.get() == NULL) {
     outer_product_.reset(
         CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
-            *A, stype, &pattern_));
+            *A, CompressedRowSparseMatrix::LOWER_TRIANGULAR, &pattern_));
   }
 
   CompressedRowSparseMatrix::ComputeOuterProduct(
-      *A, stype, pattern_, outer_product_.get());
-  cs_di lhs =
-      cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
+      *A, pattern_, outer_product_.get());
+
+  cs_di lhs = cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
 
   event_logger.AddEvent("Setup");
 
@@ -439,26 +434,22 @@
   summary.num_iterations = 1;
   summary.message = "Success.";
 
-  // Compute outerproduct to compressed row upper triangular matrix.
-  // This is the fastest option for the our default natural ordering
-  // (see comment in cholmod_factorize.c:205 in SuiteSparse).
-  const int stype = -1;
-
-  // Compute the normal equations. J'J delta = J'f and solve them
-  // using a sparse Cholesky factorization. Notice that we explicitly
-  // compute the normal equations before they can be factorized.
+  // Compute outer product to compressed row upper triangular matrix,
+  // this will be mapped to a compressed-column lower triangular
+  // matrix, which is the fastest option for our default natural
+  // ordering (see comment in cholmod_factorize.c:205 in SuiteSparse).
   if (outer_product_.get() == NULL) {
     outer_product_.reset(
         CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
-            *A, stype, &pattern_));
+            *A, CompressedRowSparseMatrix::UPPER_TRIANGULAR, &pattern_));
   }
 
   CompressedRowSparseMatrix::ComputeOuterProduct(
-      *A, stype, pattern_, outer_product_.get());
+      *A, pattern_, outer_product_.get());
 
   const int num_cols = A->num_cols();
   cholmod_sparse lhs =
-      ss_.CreateSparseMatrixTransposeView(outer_product_.get(), stype);
+      ss_.CreateSparseMatrixTransposeView(outer_product_.get());
   event_logger.AddEvent("Setup");
 
   if (options_.dynamic_sparsity) {
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index b136ce2..1a0927b 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -96,7 +96,7 @@
 }
 
 cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView(
-    CompressedRowSparseMatrix* A, const int stype) {
+    CompressedRowSparseMatrix* A) {
   cholmod_sparse m;
   m.nrow = A->num_cols();
   m.ncol = A->num_rows();
@@ -106,7 +106,15 @@
   m.i = reinterpret_cast<void*>(A->mutable_cols());
   m.x = reinterpret_cast<void*>(A->mutable_values());
   m.z = NULL;
-  m.stype = stype;
+
+  if (A->storage_type() == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+    m.stype = 1;
+  } else if (A->storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+    m.stype = -1;
+  } else {
+    m.stype = 0;
+  }
+
   m.itype = CHOLMOD_INT;
   m.xtype = CHOLMOD_REAL;
   m.dtype = CHOLMOD_DOUBLE;
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 380c16b..380d76e 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -96,11 +96,8 @@
 
   // Create a cholmod_sparse wrapper around the contents of A. This is
   // a shallow object, which refers to the contents of A and does not
-  // use the SuiteSparse machinery to allocate memory. This can
-  // create a choldmod_sparse structure with different stype (upper
-  // triangular: stype= -1 or lower triangular: stype = 1).
-  cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A,
-                                                 const int stype);
+  // use the SuiteSparse machinery to allocate memory.
+  cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A);
 
   // Given a vector x, build a cholmod_dense vector of size out_size
   // with the first in_size entries copied from x. If x is NULL, then