Cleanup reorder_program.cc
Program::SetParameterOffsetsAndIndex() was being called willy nilly.
Now the invariant is that any function that actually reorders the
program, updates the offsets and indices.
Also the logic around handling EIGEN_SPARSE has been simplified in
anticipation of the block AMD code that is forthcoming.
Last but not the least, num_eliminate_blocks, which is a rather
cryptic name to begin with has been replaced by the more meaningful
size_of_first_elimination_group.
Change-Id: I77e684f699a93b53e76aa406d64f40f8704df813
diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc
index 9b40885..4caa01f 100644
--- a/internal/ceres/reorder_program.cc
+++ b/internal/ceres/reorder_program.cc
@@ -41,7 +41,6 @@
#include "ceres/parameter_block_ordering.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
-#include "ceres/program.h"
#include "ceres/residual_block.h"
#include "ceres/solver.h"
#include "ceres/suitesparse.h"
@@ -53,12 +52,13 @@
namespace internal {
namespace {
-// Find the minimum index of any parameter block to the given residual.
-// Parameter blocks that have indices greater than num_eliminate_blocks are
-// considered to have an index equal to num_eliminate_blocks.
+// Find the minimum index of any parameter block to the given
+// residual. Parameter blocks that have indices greater than
+// size_of_first_elimination_group are considered to have an index
+// equal to size_of_first_elimination_group.
static int MinParameterBlock(const ResidualBlock* residual_block,
- int num_eliminate_blocks) {
- int min_parameter_block_position = num_eliminate_blocks;
+ int size_of_first_elimination_group) {
+ int min_parameter_block_position = size_of_first_elimination_group;
for (int i = 0; i < residual_block->NumParameterBlocks(); ++i) {
ParameterBlock* parameter_block = residual_block->parameter_blocks()[i];
if (!parameter_block->IsConstant()) {
@@ -178,30 +178,32 @@
return true;
}
-bool LexicographicallyOrderResidualBlocks(const int num_eliminate_blocks,
- Program* program,
- string* error) {
- CHECK_GE(num_eliminate_blocks, 1)
+bool LexicographicallyOrderResidualBlocks(
+ const int size_of_first_elimination_group,
+ Program* program,
+ string* error) {
+ CHECK_GE(size_of_first_elimination_group, 1)
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
// Create a histogram of the number of residuals for each E block. There is an
// extra bucket at the end to catch all non-eliminated F blocks.
- vector<int> residual_blocks_per_e_block(num_eliminate_blocks + 1);
+ vector<int> residual_blocks_per_e_block(size_of_first_elimination_group + 1);
vector<ResidualBlock*>* residual_blocks = program->mutable_residual_blocks();
vector<int> min_position_per_residual(residual_blocks->size());
for (int i = 0; i < residual_blocks->size(); ++i) {
ResidualBlock* residual_block = (*residual_blocks)[i];
- int position = MinParameterBlock(residual_block, num_eliminate_blocks);
+ int position = MinParameterBlock(residual_block,
+ size_of_first_elimination_group);
min_position_per_residual[i] = position;
- DCHECK_LE(position, num_eliminate_blocks);
+ DCHECK_LE(position, size_of_first_elimination_group);
residual_blocks_per_e_block[position]++;
}
// Run a cumulative sum on the histogram, to obtain offsets to the start of
// each histogram bucket (where each bucket is for the residuals for that
// E-block).
- vector<int> offsets(num_eliminate_blocks + 1);
+ vector<int> offsets(size_of_first_elimination_group + 1);
std::partial_sum(residual_blocks_per_e_block.begin(),
residual_blocks_per_e_block.end(),
offsets.begin());
@@ -240,7 +242,7 @@
// Sanity check #1: The difference in bucket offsets should match the
// histogram sizes.
- for (int i = 0; i < num_eliminate_blocks; ++i) {
+ for (int i = 0; i < size_of_first_elimination_group; ++i) {
CHECK_EQ(residual_blocks_per_e_block[i], offsets[i + 1] - offsets[i])
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
@@ -257,11 +259,11 @@
return true;
}
+// Pre-order the columns corresponding to the schur complement if
+// possible.
void MaybeReorderSchurComplementColumnsUsingSuiteSparse(
const ParameterBlockOrdering& parameter_block_ordering,
Program* program) {
- // Pre-order the columns corresponding to the schur complement if
- // possible.
#ifndef CERES_NO_SUITESPARSE
SuiteSparse ss;
if (!SuiteSparse::IsConstrainedApproximateMinimumDegreeOrderingAvailable()) {
@@ -283,13 +285,10 @@
// parameter_blocks.size() - 1].
MapValuesToContiguousRange(constraints.size(), &constraints[0]);
- // Set the offsets and index for CreateJacobianSparsityTranspose.
- program->SetParameterOffsetsAndIndex();
// Compute a block sparse presentation of J'.
scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
program->CreateJacobianBlockSparsityTranspose());
-
cholmod_sparse* block_jacobian_transpose =
ss.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
@@ -303,6 +302,8 @@
for (int i = 0; i < program->NumParameterBlocks(); ++i) {
parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
}
+
+ program->SetParameterOffsetsAndIndex();
#endif
}
@@ -313,7 +314,8 @@
ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* error) {
- if (parameter_block_ordering->NumElements() != program->NumParameterBlocks()) {
+ if (parameter_block_ordering->NumElements() !=
+ program->NumParameterBlocks()) {
*error = StringPrintf(
"The program has %d parameter blocks, but the parameter block "
"ordering has %d parameter blocks.",
@@ -330,7 +332,7 @@
// this means that the user wishes for Ceres to identify the
// e_blocks, which we do by computing a maximal independent set.
vector<ParameterBlock*> schur_ordering;
- const int num_eliminate_blocks =
+ const int size_of_first_elimination_group =
ComputeStableSchurOrdering(*program, &schur_ordering);
CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks())
@@ -340,7 +342,7 @@
// Update the parameter_block_ordering object.
for (int i = 0; i < schur_ordering.size(); ++i) {
double* parameter_block = schur_ordering[i]->mutable_user_state();
- const int group_id = (i < num_eliminate_blocks) ? 0 : 1;
+ const int group_id = (i < size_of_first_elimination_group) ? 0 : 1;
parameter_block_ordering->AddElementToGroup(parameter_block, group_id);
}
@@ -373,6 +375,8 @@
}
}
+ program->SetParameterOffsetsAndIndex();
+
if (linear_solver_type == SPARSE_SCHUR &&
sparse_linear_algebra_library_type == SUITE_SPARSE) {
MaybeReorderSchurComplementColumnsUsingSuiteSparse(
@@ -380,18 +384,18 @@
program);
}
- program->SetParameterOffsetsAndIndex();
+
+
// Schur type solvers also require that their residual blocks be
// lexicographically ordered.
- const int num_eliminate_blocks =
+ const int size_of_first_elimination_group =
parameter_block_ordering->group_to_elements().begin()->second.size();
- if (!LexicographicallyOrderResidualBlocks(num_eliminate_blocks,
+ if (!LexicographicallyOrderResidualBlocks(size_of_first_elimination_group,
program,
error)) {
return false;
}
- program->SetParameterOffsetsAndIndex();
return true;
}
@@ -400,30 +404,6 @@
const ParameterBlockOrdering& parameter_block_ordering,
Program* program,
string* error) {
-
- if (sparse_linear_algebra_library_type != SUITE_SPARSE &&
- sparse_linear_algebra_library_type != CX_SPARSE &&
- sparse_linear_algebra_library_type != EIGEN_SPARSE) {
- *error = "Unknown sparse linear algebra library.";
- return false;
- }
-
- // For Eigen, there is nothing to do. This is because Eigen in its
- // current stable version does not expose a method for doing
- // symbolic analysis on pre-ordered matrices, so a block
- // pre-ordering is a bit pointless.
- //
- // The dev version as recently as July 20, 2014 has support for
- // pre-ordering. Once this becomes more widespread, or we add
- // support for detecting Eigen versions, we can add support for this
- // along the lines of CXSparse.
- if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
- program->SetParameterOffsetsAndIndex();
- return true;
- }
-
- // Set the offsets and index for CreateJacobianSparsityTranspose.
- program->SetParameterOffsetsAndIndex();
// Compute a block sparse presentation of J'.
scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
program->CreateJacobianBlockSparsityTranspose());
@@ -438,10 +418,16 @@
parameter_blocks,
parameter_block_ordering,
&ordering[0]);
- } else if (sparse_linear_algebra_library_type == CX_SPARSE){
+ } else if (sparse_linear_algebra_library_type == CX_SPARSE) {
OrderingForSparseNormalCholeskyUsingCXSparse(
*tsm_block_jacobian_transpose,
&ordering[0]);
+ } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
+ // Starting with v3.2.2 Eigen has support for symbolic analysis on
+ // pre-ordered matrices.
+ //
+ // TODO(sameeragarwal): Apply block amd for eigen.
+ return true;
}
// Apply ordering.
diff --git a/internal/ceres/reorder_program.h b/internal/ceres/reorder_program.h
index d3962f9..0474c1f 100644
--- a/internal/ceres/reorder_program.h
+++ b/internal/ceres/reorder_program.h
@@ -51,7 +51,7 @@
// Reorder the residuals for program, if necessary, so that the residuals
// involving each E block occur together. This is a necessary condition for the
// Schur eliminator, which works on these "row blocks" in the jacobian.
-bool LexicographicallyOrderResidualBlocks(int num_eliminate_blocks,
+bool LexicographicallyOrderResidualBlocks(int size_of_first_elimination_group,
Program* program,
string* error);
