Block ordering for SPARSE_SCHUR + CX_SPARSE. Uptil now only SuiteSparse when used with SPARSE_SCHUR would use the block structure of the reduced camera matrix to find a fill-reducing ordering. This leads to substantial speedup for some bundle adjustment problems. Credit for this technique goes to the authors of g2o. I learned about it from reading their source code. Change-Id: I5403efefd4d9552c9c6fc6e02a65498bdf171584
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 3734c08..549c94e 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -40,6 +40,7 @@ block_structure.cc canonical_views_clustering.cc cgnr_solver.cc + compressed_col_sparse_matrix_utils.cc compressed_row_jacobian_writer.cc compressed_row_sparse_matrix.cc conditioned_cost_function.cc @@ -280,8 +281,10 @@ CERES_TEST(schur_eliminator) CERES_TEST(solver_impl) + # TODO(sameeragarwal): This test should ultimately be made + # independent of SuiteSparse. IF (${SUITESPARSE_FOUND}) - CERES_TEST(suitesparse) + CERES_TEST(compressed_col_sparse_matrix_utils) ENDIF (${SUITESPARSE_FOUND}) CERES_TEST(symmetric_linear_solver)
diff --git a/internal/ceres/compressed_col_sparse_matrix_utils.cc b/internal/ceres/compressed_col_sparse_matrix_utils.cc new file mode 100644 index 0000000..ba76dad --- /dev/null +++ b/internal/ceres/compressed_col_sparse_matrix_utils.cc
@@ -0,0 +1,117 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2013 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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/compressed_col_sparse_matrix_utils.h" + +#include <vector> +#include "ceres/internal/port.h" +#include "glog/logging.h" + +namespace ceres { +namespace internal { + +void CompressedColumnScalarMatrixToBlockMatrix(const int* scalar_rows, + const int* scalar_cols, + 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]; + } + + // This loop extracts the block sparsity of the scalar sparse matrix + // 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 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++; + } + } +} +} // namespace internal +} // namespace ceres
diff --git a/internal/ceres/compressed_col_sparse_matrix_utils.h b/internal/ceres/compressed_col_sparse_matrix_utils.h new file mode 100644 index 0000000..afabf1c --- /dev/null +++ b/internal/ceres/compressed_col_sparse_matrix_utils.h
@@ -0,0 +1,67 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2013 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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_COMPRESSED_COL_SPARSE_MATRIX_UTILS_H_ +#define CERES_INTERNAL_COMPRESSED_COL_SPARSE_MATRIX_UTILS_H_ + +#include <vector> +#include "ceres/internal/port.h" + +namespace ceres { +namespace internal { + +// Extract the block sparsity pattern of the scalar compressed columns +// matrix and return it in compressed column form. The compressed +// column form is stored in two vectors block_rows, and block_cols, +// which correspond to the row and column arrays in a compressed +// column sparse matrix. +// +// If c_ij is the block in the matrix A corresponding to row block i +// and column block j, then it is expected that A contains at least +// one non-zero entry corresponding to the top left entry of c_ij, +// as that entry is used to detect the presence of a non-zero c_ij. +void CompressedColumnScalarMatrixToBlockMatrix(const int* scalar_rows, + const int* scalar_cols, + const vector<int>& row_blocks, + const vector<int>& col_blocks, + vector<int>* block_rows, + vector<int>* block_cols); + +// Given a set of blocks and a permutation of these blocks, compute +// the corresponding "scalar" ordering, where the scalar ordering of +// size sum(blocks). +void BlockOrderingToScalarOrdering(const vector<int>& blocks, + const vector<int>& block_ordering, + vector<int>* scalar_ordering); + +} // namespace internal +} // namespace ceres + +#endif // CERES_INTERNAL_COMPRESSED_COL_SPARSE_MATRIX_UTILS_H_
diff --git a/internal/ceres/suitesparse_test.cc b/internal/ceres/compressed_col_sparse_matrix_utils_test.cc similarity index 88% rename from internal/ceres/suitesparse_test.cc rename to internal/ceres/compressed_col_sparse_matrix_utils_test.cc index 72e3e68..7efa0e3 100644 --- a/internal/ceres/suitesparse_test.cc +++ b/internal/ceres/compressed_col_sparse_matrix_utils_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2012 Google Inc. All rights reserved. +// Copyright 2013 Google Inc. All rights reserved. // http://code.google.com/p/ceres-solver/ // // Redistribution and use in source and binary forms, with or without @@ -29,6 +29,7 @@ // Author: sameeragarwal@google.com (Sameer Agarwal) #include <algorithm> +#include "ceres/compressed_col_sparse_matrix_utils.h" #include "ceres/internal/port.h" #include "ceres/suitesparse.h" #include "ceres/triplet_sparse_matrix.h" @@ -38,7 +39,7 @@ namespace ceres { namespace internal { -TEST(SuiteSparse, BlockPermutationToScalarPermutation) { +TEST(_, BlockPermutationToScalarPermutation) { vector<int> blocks; // Block structure // 0 --1- ---2--- ---3--- 4 @@ -73,9 +74,9 @@ expected_scalar_ordering.push_back(8); vector<int> scalar_ordering; - SuiteSparse::BlockOrderingToScalarOrdering(blocks, - block_ordering, - &scalar_ordering); + BlockOrderingToScalarOrdering(blocks, + block_ordering, + &scalar_ordering); EXPECT_EQ(scalar_ordering.size(), expected_scalar_ordering.size()); for (int i = 0; i < expected_scalar_ordering.size(); ++i) { EXPECT_EQ(scalar_ordering[i], expected_scalar_ordering[i]); @@ -109,7 +110,7 @@ return offset; } -TEST(SuiteSparse, ScalarMatrixToBlockMatrix) { +TEST(_, ScalarMatrixToBlockMatrix) { // Block sparsity. // // [1 2 3 2] @@ -170,11 +171,12 @@ vector<int> block_rows; vector<int> block_cols; - SuiteSparse::ScalarMatrixToBlockMatrix(ccsm.get(), - row_blocks, - col_blocks, - &block_rows, - &block_cols); + CompressedColumnScalarMatrixToBlockMatrix(reinterpret_cast<const int*>(ccsm->i), + reinterpret_cast<const int*>(ccsm->p), + row_blocks, + col_blocks, + &block_rows, + &block_cols); EXPECT_EQ(block_cols.size(), expected_block_cols.size()); EXPECT_EQ(block_rows.size(), expected_block_rows.size());
diff --git a/internal/ceres/cxsparse.cc b/internal/ceres/cxsparse.cc index 3fbc271..b7f2520 100644 --- a/internal/ceres/cxsparse.cc +++ b/internal/ceres/cxsparse.cc
@@ -32,7 +32,10 @@ #include "ceres/cxsparse.h" +#include <vector> +#include "ceres/compressed_col_sparse_matrix_utils.h" #include "ceres/compressed_row_sparse_matrix.h" +#include "ceres/internal/port.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" @@ -44,24 +47,25 @@ CXSparse::~CXSparse() { if (scratch_size_ > 0) { - cs_free(scratch_); + cs_di_free(scratch_); } } + bool CXSparse::SolveCholesky(cs_di* A, cs_dis* symbolic_factorization, double* b) { // Make sure we have enough scratch space available. if (scratch_size_ < A->n) { if (scratch_size_ > 0) { - cs_free(scratch_); + cs_di_free(scratch_); } - scratch_ = reinterpret_cast<CS_ENTRY*>(cs_malloc(A->n, sizeof(CS_ENTRY))); + scratch_ = reinterpret_cast<CS_ENTRY*>(cs_di_malloc(A->n, sizeof(CS_ENTRY))); scratch_size_ = A->n; } // Solve using Cholesky factorization - csn* numeric_factorization = cs_chol(A, symbolic_factorization); + csn* numeric_factorization = cs_di_chol(A, symbolic_factorization); if (numeric_factorization == NULL) { LOG(WARNING) << "Cholesky factorization failed."; return false; @@ -71,19 +75,16 @@ // succeeded as well. In the comments below, "x" refers to the scratch space. // // Set x = P * b. - cs_ipvec(symbolic_factorization->pinv, b, scratch_, A->n); - + cs_di_ipvec(symbolic_factorization->pinv, b, scratch_, A->n); // Set x = L \ x. - cs_lsolve(numeric_factorization->L, scratch_); - + cs_di_lsolve(numeric_factorization->L, scratch_); // Set x = L' \ x. - cs_ltsolve(numeric_factorization->L, scratch_); - + cs_di_ltsolve(numeric_factorization->L, scratch_); // Set b = P' * x. - cs_pvec(symbolic_factorization->pinv, scratch_, b, A->n); + cs_di_pvec(symbolic_factorization->pinv, scratch_, b, A->n); // Free Cholesky factorization. - cs_nfree(numeric_factorization); + cs_di_nfree(numeric_factorization); return true; } @@ -92,6 +93,61 @@ return cs_schol(1, A); } +cs_dis* CXSparse::BlockAnalyzeCholesky(cs_di* A, + const vector<int>& row_blocks, + const vector<int>& col_blocks) { + const int num_row_blocks = row_blocks.size(); + const int num_col_blocks = col_blocks.size(); + + vector<int> block_rows; + vector<int> block_cols; + CompressedColumnScalarMatrixToBlockMatrix(A->i, + A->p, + row_blocks, + col_blocks, + &block_rows, + &block_cols); + cs_di block_matrix; + block_matrix.m = num_row_blocks; + block_matrix.n = num_col_blocks; + block_matrix.nz = -1; + block_matrix.nzmax = block_rows.size(); + block_matrix.p = &block_cols[0]; + block_matrix.i = &block_rows[0]; + block_matrix.x = NULL; + + int* ordering = cs_amd(1, &block_matrix); + vector<int> block_ordering(num_row_blocks, -1); + copy(ordering, ordering + num_row_blocks, &block_ordering[0]); + cs_free(ordering); + + vector<int> scalar_ordering; + BlockOrderingToScalarOrdering(row_blocks, block_ordering, &scalar_ordering); + + cs_dis* symbolic_factorization = reinterpret_cast<cs_dis*>(cs_calloc(1, sizeof(cs_dis))); + symbolic_factorization->pinv = cs_pinv(&scalar_ordering[0], A->n); + cs* permuted_A = cs_symperm(A, symbolic_factorization->pinv, 0); + + symbolic_factorization->parent = cs_etree(permuted_A, 0); + int* postordering = cs_post(symbolic_factorization->parent, A->n); + int* column_counts = cs_counts(permuted_A, symbolic_factorization->parent, postordering, 0); + cs_free(postordering); + cs_spfree(permuted_A); + + symbolic_factorization->cp = (int*) cs_malloc(A->n+1, sizeof(int)); + symbolic_factorization->lnz = cs_cumsum(symbolic_factorization->cp, column_counts, A->n); + symbolic_factorization->unz = symbolic_factorization->lnz; + + cs_free(column_counts); + + if (symbolic_factorization->lnz < 0) { + cs_sfree(symbolic_factorization); + symbolic_factorization = NULL; + } + + return symbolic_factorization; +} + cs_di CXSparse::CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A) { cs_di At; At.m = A->num_cols();
diff --git a/internal/ceres/cxsparse.h b/internal/ceres/cxsparse.h index dd5eadc..d34b635 100644 --- a/internal/ceres/cxsparse.h +++ b/internal/ceres/cxsparse.h
@@ -33,7 +33,9 @@ #ifndef CERES_NO_CXSPARSE +#include <vector> #include "cs.h" +#include "ceres/internal/port.h" namespace ceres { namespace internal { @@ -72,6 +74,10 @@ // The returned matrix should be deallocated with Free when not used anymore. cs_dis* AnalyzeCholesky(cs_di* A); + cs_dis* BlockAnalyzeCholesky(cs_di* A, + const vector<int>& row_blocks, + const vector<int>& col_blocks); + void Free(cs_di* sparse_matrix); void Free(cs_dis* symbolic_factorization);
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index 9539c0d..0defcd6 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc
@@ -339,7 +339,7 @@ // Compute symbolic factorization if not available. if (cxsparse_factor_ == NULL) { - cxsparse_factor_ = CHECK_NOTNULL(cxsparse_.AnalyzeCholesky(lhs)); + cxsparse_factor_ = CHECK_NOTNULL(cxsparse_.BlockAnalyzeCholesky(lhs, blocks_, blocks_)); } // Solve the linear system.
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc index d576f6a..fe2edd3 100644 --- a/internal/ceres/suitesparse.cc +++ b/internal/ceres/suitesparse.cc
@@ -33,6 +33,7 @@ #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" @@ -202,11 +203,12 @@ 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; @@ -231,88 +233,6 @@ 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);
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h index a1a4f35..e138623 100644 --- a/internal/ceres/suitesparse.h +++ b/internal/ceres/suitesparse.h
@@ -184,29 +184,6 @@ const vector<int>& col_blocks, vector<int>* ordering); - // Given a set of blocks and a permutation of these blocks, compute - // the corresponding "scalar" ordering, where the scalar ordering of - // size sum(blocks). - static void BlockOrderingToScalarOrdering(const vector<int>& blocks, - const vector<int>& block_ordering, - vector<int>* scalar_ordering); - - // Extract the block sparsity pattern of the scalar sparse matrix - // A and return it in compressed column form. The compressed column - // form is stored in two vectors block_rows, and block_cols, which - // correspond to the row and column arrays in a compressed column sparse - // matrix. - // - // If c_ij is the block in the matrix A corresponding to row block i - // and column block j, then it is expected that A contains at least - // one non-zero entry corresponding to the top left entry of c_ij, - // as that entry is used to detect the presence of a non-zero c_ij. - static void ScalarMatrixToBlockMatrix(const cholmod_sparse* A, - const vector<int>& row_blocks, - const vector<int>& col_blocks, - vector<int>* block_rows, - vector<int>* block_cols); - // Find a fill reducing approximate minimum degree // ordering. ordering is expected to be large enough to hold the // ordering.