Add a storage type to CompressedRowSparseMatrix

By adding an enum to CompressedRowSparseMatrix, which indicates
whether the matrix is unsymmetric, upper or lower triangular
we are able to improve the readability and fix some minor
bugs in the way some matrix manipulation code was being
called.

Thank to William Rucklidge for this suggestion.

Change-Id: I355c90d11cd5d31f5a25741b0bda4fc4583e9095

internal/ceres/compressed_row_sparse_matrix.cc

@@ -109,6 +109,7 @@ CompressedRowSparseMatrix::CompressedRowSparseMatrix(int num_rows,
                                                      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 @@ CompressedRowSparseMatrix::CompressedRowSparseMatrix(
     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 @@ CompressedRowSparseMatrix::CompressedRowSparseMatrix(const double* diagonal,
 
   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* CompressedRowSparseMatrix::Transpose() const {
   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 @@ struct ProductTerm {
 // 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* CreateOuterProductMatrix(
 
   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* CreateOuterProductMatrix(
 
 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 @@ CompressedRowSparseMatrix* CompressAndFillProgram(
 
   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 @@ void ComputeBlockMultiplication(const int row_block_size,
 
 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 @@ CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
   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 @@ CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
   }
 
   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 @@ CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
 // 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 @@ void CompressedRowSparseMatrix::ComputeOuterProduct(
   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 @@ void CompressedRowSparseMatrix::ComputeOuterProduct(
       // 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 @@ void CompressedRowSparseMatrix::ComputeOuterProduct(
                                      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) {

internal/ceres/compressed_row_sparse_matrix.h

@@ -48,6 +48,13 @@ class TripletSparseMatrix;
 
 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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
 
   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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   // 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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   // 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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   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 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   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

internal/ceres/compressed_row_sparse_matrix_test.cc

@@ -494,16 +494,16 @@ TEST(CompressedRowSparseMatrix, ComputeOuterProduct) {
         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>(

internal/ceres/sparse_normal_cholesky_solver.cc

@@ -275,21 +275,18 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingEigen(
                                &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 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
   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 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
   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) {

internal/ceres/suitesparse.cc

@@ -96,7 +96,7 @@ cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose(
 }
 
 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 @@ cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView(
   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;

internal/ceres/suitesparse.h

@@ -96,11 +96,8 @@ class SuiteSparse {
 
   // 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