diff options
Diffstat (limited to 'internal/ceres/suitesparse.cc')
| -rw-r--r-- | internal/ceres/suitesparse.cc | 224 |
1 files changed, 111 insertions, 113 deletions
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc index cf3c48f..9de32fd 100644 --- a/internal/ceres/suitesparse.cc +++ b/internal/ceres/suitesparse.cc @@ -33,10 +33,21 @@ #include <vector> #include "cholmod.h" +#include "ceres/compressed_col_sparse_matrix_utils.h" #include "ceres/compressed_row_sparse_matrix.h" #include "ceres/triplet_sparse_matrix.h" + namespace ceres { namespace internal { + +SuiteSparse::SuiteSparse() { + cholmod_start(&cc_); +} + +SuiteSparse::~SuiteSparse() { + cholmod_finish(&cc_); +} + cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) { cholmod_triplet triplet; @@ -77,23 +88,23 @@ cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose( return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_); } -cholmod_sparse* SuiteSparse::CreateSparseMatrixTransposeView( +cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView( CompressedRowSparseMatrix* A) { - cholmod_sparse* m = new cholmod_sparse_struct; - m->nrow = A->num_cols(); - m->ncol = A->num_rows(); - m->nzmax = A->num_nonzeros(); - - m->p = reinterpret_cast<void*>(A->mutable_rows()); - m->i = reinterpret_cast<void*>(A->mutable_cols()); - m->x = reinterpret_cast<void*>(A->mutable_values()); - - m->stype = 0; // Matrix is not symmetric. - m->itype = CHOLMOD_INT; - m->xtype = CHOLMOD_REAL; - m->dtype = CHOLMOD_DOUBLE; - m->sorted = 1; - m->packed = 1; + cholmod_sparse m; + m.nrow = A->num_cols(); + m.ncol = A->num_rows(); + m.nzmax = A->num_nonzeros(); + m.nz = NULL; + m.p = reinterpret_cast<void*>(A->mutable_rows()); + m.i = reinterpret_cast<void*>(A->mutable_cols()); + m.x = reinterpret_cast<void*>(A->mutable_values()); + m.z = NULL; + m.stype = 0; // Matrix is not symmetric. + m.itype = CHOLMOD_INT; + m.xtype = CHOLMOD_REAL; + m.dtype = CHOLMOD_DOUBLE; + m.sorted = 1; + m.packed = 1; return m; } @@ -117,10 +128,16 @@ cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A) { cc_.nmethods = 1; cc_.method[0].ordering = CHOLMOD_AMD; cc_.supernodal = CHOLMOD_AUTO; + cholmod_factor* factor = cholmod_analyze(A, &cc_); CHECK_EQ(cc_.status, CHOLMOD_OK) << "Cholmod symbolic analysis failed " << cc_.status; CHECK_NOTNULL(factor); + + if (VLOG_IS_ON(2)) { + cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); + } + return factor; } @@ -135,16 +152,42 @@ cholmod_factor* SuiteSparse::BlockAnalyzeCholesky( return AnalyzeCholeskyWithUserOrdering(A, ordering); } -cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(cholmod_sparse* A, - const vector<int>& ordering) { +cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering( + cholmod_sparse* A, + const vector<int>& ordering) { CHECK_EQ(ordering.size(), A->nrow); - cc_.nmethods = 1 ; + + cc_.nmethods = 1; cc_.method[0].ordering = CHOLMOD_GIVEN; + cholmod_factor* factor = cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), NULL, 0, &cc_); CHECK_EQ(cc_.status, CHOLMOD_OK) << "Cholmod symbolic analysis failed " << cc_.status; CHECK_NOTNULL(factor); + + if (VLOG_IS_ON(2)) { + cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); + } + + return factor; +} + +cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering( + cholmod_sparse* A) { + cc_.nmethods = 1; + cc_.method[0].ordering = CHOLMOD_NATURAL; + cc_.postorder = 0; + + cholmod_factor* factor = cholmod_analyze(A, &cc_); + CHECK_EQ(cc_.status, CHOLMOD_OK) + << "Cholmod symbolic analysis failed " << cc_.status; + CHECK_NOTNULL(factor); + + if (VLOG_IS_ON(2)) { + cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); + } + return factor; } @@ -160,11 +203,12 @@ bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A, vector<int> block_cols; vector<int> block_rows; - ScalarMatrixToBlockMatrix(A, - row_blocks, - col_blocks, - &block_rows, - &block_cols); + CompressedColumnScalarMatrixToBlockMatrix(reinterpret_cast<const int*>(A->i), + reinterpret_cast<const int*>(A->p), + row_blocks, + col_blocks, + &block_rows, + &block_cols); cholmod_sparse_struct block_matrix; block_matrix.nrow = num_row_blocks; @@ -189,122 +233,56 @@ bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A, return true; } -void SuiteSparse::ScalarMatrixToBlockMatrix(const cholmod_sparse* A, - const vector<int>& row_blocks, - const vector<int>& col_blocks, - vector<int>* block_rows, - vector<int>* block_cols) { - CHECK_NOTNULL(block_rows)->clear(); - CHECK_NOTNULL(block_cols)->clear(); - const int num_row_blocks = row_blocks.size(); - const int num_col_blocks = col_blocks.size(); - - vector<int> row_block_starts(num_row_blocks); - for (int i = 0, cursor = 0; i < num_row_blocks; ++i) { - row_block_starts[i] = cursor; - cursor += row_blocks[i]; - } - - // The reinterpret_cast is needed here because CHOLMOD stores arrays - // as void*. - const int* scalar_cols = reinterpret_cast<const int*>(A->p); - const int* scalar_rows = reinterpret_cast<const int*>(A->i); - - // This loop extracts the block sparsity of the scalar sparse matrix - // A. It does so by iterating over the columns, but only considering - // the columns corresponding to the first element of each column - // block. Within each column, the inner loop iterates over the rows, - // and detects the presence of a row block by checking for the - // presence of a non-zero entry corresponding to its first element. - block_cols->push_back(0); - int c = 0; - for (int col_block = 0; col_block < num_col_blocks; ++col_block) { - int column_size = 0; - for (int idx = scalar_cols[c]; idx < scalar_cols[c + 1]; ++idx) { - vector<int>::const_iterator it = lower_bound(row_block_starts.begin(), - row_block_starts.end(), - scalar_rows[idx]); - // Since we are using lower_bound, it will return the row id - // where the row block starts. For everything but the first row - // of the block, where these values will be the same, we can - // skip, as we only need the first row to detect the presence of - // the block. - // - // For rows all but the first row in the last row block, - // lower_bound will return row_block_starts.end(), but those can - // be skipped like the rows in other row blocks too. - if (it == row_block_starts.end() || *it != scalar_rows[idx]) { - continue; - } - - block_rows->push_back(it - row_block_starts.begin()); - ++column_size; - } - block_cols->push_back(block_cols->back() + column_size); - c += col_blocks[col_block]; - } -} - -void SuiteSparse::BlockOrderingToScalarOrdering( - const vector<int>& blocks, - const vector<int>& block_ordering, - vector<int>* scalar_ordering) { - CHECK_EQ(blocks.size(), block_ordering.size()); - const int num_blocks = blocks.size(); - - // block_starts = [0, block1, block1 + block2 ..] - vector<int> block_starts(num_blocks); - for (int i = 0, cursor = 0; i < num_blocks ; ++i) { - block_starts[i] = cursor; - cursor += blocks[i]; - } - - scalar_ordering->resize(block_starts.back() + blocks.back()); - int cursor = 0; - for (int i = 0; i < num_blocks; ++i) { - const int block_id = block_ordering[i]; - const int block_size = blocks[block_id]; - int block_position = block_starts[block_id]; - for (int j = 0; j < block_size; ++j) { - (*scalar_ordering)[cursor++] = block_position++; - } - } -} - bool SuiteSparse::Cholesky(cholmod_sparse* A, cholmod_factor* L) { CHECK_NOTNULL(A); CHECK_NOTNULL(L); + // Save the current print level and silence CHOLMOD, otherwise + // CHOLMOD is prone to dumping stuff to stderr, which can be + // distracting when the error (matrix is indefinite) is not a fatal + // failure. + const int old_print_level = cc_.print; + cc_.print = 0; + cc_.quick_return_if_not_posdef = 1; int status = cholmod_factorize(A, L, &cc_); + cc_.print = old_print_level; + + // TODO(sameeragarwal): This switch statement is not consistent. It + // treats all kinds of CHOLMOD failures as warnings. Some of these + // like out of memory are definitely not warnings. The problem is + // that the return value Cholesky is two valued, but the state of + // the linear solver is really three valued. SUCCESS, + // NON_FATAL_FAILURE (e.g., indefinite matrix) and FATAL_FAILURE + // (e.g. out of memory). switch (cc_.status) { case CHOLMOD_NOT_INSTALLED: - LOG(WARNING) << "Cholmod failure: method not installed."; + LOG(WARNING) << "CHOLMOD failure: Method not installed."; return false; case CHOLMOD_OUT_OF_MEMORY: - LOG(WARNING) << "Cholmod failure: out of memory."; + LOG(WARNING) << "CHOLMOD failure: Out of memory."; return false; case CHOLMOD_TOO_LARGE: - LOG(WARNING) << "Cholmod failure: integer overflow occured."; + LOG(WARNING) << "CHOLMOD failure: Integer overflow occured."; return false; case CHOLMOD_INVALID: - LOG(WARNING) << "Cholmod failure: invalid input."; + LOG(WARNING) << "CHOLMOD failure: Invalid input."; return false; case CHOLMOD_NOT_POSDEF: // TODO(sameeragarwal): These two warnings require more // sophisticated handling going forward. For now we will be // strict and treat them as failures. - LOG(WARNING) << "Cholmod warning: matrix not positive definite."; + LOG(WARNING) << "CHOLMOD warning: Matrix not positive definite."; return false; case CHOLMOD_DSMALL: - LOG(WARNING) << "Cholmod warning: D for LDL' or diag(L) or " + LOG(WARNING) << "CHOLMOD warning: D for LDL' or diag(L) or " << "LL' has tiny absolute value."; return false; case CHOLMOD_OK: if (status != 0) { return true; } - LOG(WARNING) << "Cholmod failure: cholmod_factorize returned zero " + LOG(WARNING) << "CHOLMOD failure: cholmod_factorize returned zero " << "but cholmod_common::status is CHOLMOD_OK." << "Please report this to ceres-solver@googlegroups.com."; return false; @@ -340,6 +318,26 @@ cholmod_dense* SuiteSparse::SolveCholesky(cholmod_sparse* A, return NULL; } +void SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, + int* ordering) { + cholmod_amd(matrix, NULL, 0, ordering, &cc_); +} + +void SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering( + cholmod_sparse* matrix, + int* constraints, + int* ordering) { +#ifndef CERES_NO_CAMD + cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_); +#else + LOG(FATAL) << "Congratulations you have found a bug in Ceres." + << "Ceres Solver was compiled with SuiteSparse " + << "version 4.1.0 or less. Calling this function " + << "in that case is a bug. Please contact the" + << "the Ceres Solver developers."; +#endif +} + } // namespace internal } // namespace ceres |