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);