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.
