Fix a bug in the Schur eliminator
The schur eliminator treats rows with e blocks and row with
no e blocks separately. The template specialization logic only
applies to the rows with e blocks.
So, in cases where the rows with e-blocks have a fixed size f-block
but the rows without e-blocks have f-blocks of varying sizes,
DetectStructure will return a static f-block size, but we need to be
careful that we do not blindly use that static f-block size everywhere.
This patch fixes a bug where such care was not being taken, where
it was assumed that the static f-block size could be assumed for all
f-block sizes.
A new test is added, which triggers an exception in debug mode. In
release mode this error does not present itself, due to a peculiarity
of the way Eigen works.
Thanks to Werner Trobin for reporting this bug.
Change-Id: I8ae7aabf8eed8c3f9cf74b6c74d632ba44f82581
diff --git a/internal/ceres/schur_eliminator_impl.h b/internal/ceres/schur_eliminator_impl.h
index 7e06806..f253588 100644
--- a/internal/ceres/schur_eliminator_impl.h
+++ b/internal/ceres/schur_eliminator_impl.h
@@ -194,7 +194,7 @@
&row_stride, &col_stride);
if (cell_info != NULL) {
const int block_size = bs->cols[i].size;
- typename EigenTypes<kFBlockSize>::ConstVectorRef
+ typename EigenTypes<Eigen::Dynamic>::ConstVectorRef
diag(D + bs->cols[i].position, block_size);
CeresMutexLock l(&cell_info->m);
@@ -205,20 +205,19 @@
}
}
- // Eliminate y blocks one chunk at a time. For each chunk,x3
- // compute the entries of the normal equations and the gradient
- // vector block corresponding to the y block and then apply
- // Gaussian elimination to them. The matrix ete stores the normal
- // matrix corresponding to the block being eliminated and array
- // buffer_ contains the non-zero blocks in the row corresponding
- // to this y block in the normal equations. This computation is
- // done in ChunkDiagonalBlockAndGradient. UpdateRhs then applies
- // gaussian elimination to the rhs of the normal equations,
- // updating the rhs of the reduced linear system by modifying rhs
- // blocks for all the z blocks that share a row block/residual
- // term with the y block. EliminateRowOuterProduct does the
- // corresponding operation for the lhs of the reduced linear
- // system.
+ // Eliminate y blocks one chunk at a time. For each chunk, compute
+ // the entries of the normal equations and the gradient vector block
+ // corresponding to the y block and then apply Gaussian elimination
+ // to them. The matrix ete stores the normal matrix corresponding to
+ // the block being eliminated and array buffer_ contains the
+ // non-zero blocks in the row corresponding to this y block in the
+ // normal equations. This computation is done in
+ // ChunkDiagonalBlockAndGradient. UpdateRhs then applies gaussian
+ // elimination to the rhs of the normal equations, updating the rhs
+ // of the reduced linear system by modifying rhs blocks for all the
+ // z blocks that share a row block/residual term with the y
+ // block. EliminateRowOuterProduct does the corresponding operation
+ // for the lhs of the reduced linear system.
#pragma omp parallel for num_threads(num_threads_) schedule(dynamic)
for (int i = 0; i < chunks_.size(); ++i) {
#ifdef CERES_USE_OPENMP
@@ -594,8 +593,8 @@
void
SchurEliminator<kRowBlockSize, kEBlockSize, kFBlockSize>::
NoEBlockRowOuterProduct(const BlockSparseMatrix* A,
- int row_block_index,
- BlockRandomAccessMatrix* lhs) {
+ int row_block_index,
+ BlockRandomAccessMatrix* lhs) {
const CompressedRowBlockStructure* bs = A->block_structure();
const CompressedRow& row = bs->rows[row_block_index];
const double* values = A->values();
