Relax the limitation that SchurEliminator::Eliminate requires a rhs.
When using the SchurEliminator to compute a preconditioner, there
is no rhs. This CL removes that limitation and simplifies the
call sites in the two preconditioners.
https://github.com/ceres-solver/ceres-solver/issues/271
Change-Id: I05b8518fdc9f0d1a6d88ae76d3a7e8e838e7204a
diff --git a/internal/ceres/schur_eliminator.h b/internal/ceres/schur_eliminator.h
index d93bcb0..deae5a5 100644
--- a/internal/ceres/schur_eliminator.h
+++ b/internal/ceres/schur_eliminator.h
@@ -218,9 +218,6 @@
// parameters as template arguments so that they are visible to the
// compiler and can be used for compile time optimization of the low
// level linear algebra routines.
-//
-// This implementation is mulithreaded using OpenMP. The level of
-// parallelism is controlled by LinearSolver::Options::num_threads.
template <int kRowBlockSize = Eigen::Dynamic,
int kEBlockSize = Eigen::Dynamic,
int kFBlockSize = Eigen::Dynamic >
diff --git a/internal/ceres/schur_eliminator_impl.h b/internal/ceres/schur_eliminator_impl.h
index 06d0983..bc74015 100644
--- a/internal/ceres/schur_eliminator_impl.h
+++ b/internal/ceres/schur_eliminator_impl.h
@@ -176,7 +176,9 @@
double* rhs) {
if (lhs->num_rows() > 0) {
lhs->SetZero();
- VectorRef(rhs, lhs->num_rows()).setZero();
+ if (rhs) {
+ VectorRef(rhs, lhs->num_rows()).setZero();
+ }
}
const CompressedRowBlockStructure* bs = A->block_structure();
@@ -276,15 +278,16 @@
//
// rhs = F'b - F'E(E'E)^(-1) E'b
- FixedArray<double, 8> inverse_ete_g(e_block_size);
- MatrixVectorMultiply<kEBlockSize, kEBlockSize, 0>(
- inverse_ete.data(),
- e_block_size,
- e_block_size,
- g.get(),
- inverse_ete_g.get());
-
- UpdateRhs(chunk, A, b, chunk.start, inverse_ete_g.get(), rhs);
+ if (rhs) {
+ FixedArray<double, 8> inverse_ete_g(e_block_size);
+ MatrixVectorMultiply<kEBlockSize, kEBlockSize, 0>(
+ inverse_ete.data(),
+ e_block_size,
+ e_block_size,
+ g.get(),
+ inverse_ete_g.get());
+ UpdateRhs(chunk, A, b, chunk.start, inverse_ete_g.get(), rhs);
+ }
// S -= F'E(E'E)^{-1}E'F
ChunkOuterProduct(
@@ -469,12 +472,13 @@
values + e_cell.position, row.block.size, e_block_size,
ete->data(), 0, 0, e_block_size, e_block_size);
- // g += E_i' b_i
- MatrixTransposeVectorMultiply<kRowBlockSize, kEBlockSize, 1>(
- values + e_cell.position, row.block.size, e_block_size,
- b + b_pos,
- g);
-
+ if (b) {
+ // g += E_i' b_i
+ MatrixTransposeVectorMultiply<kRowBlockSize, kEBlockSize, 1>(
+ values + e_cell.position, row.block.size, e_block_size,
+ b + b_pos,
+ g);
+ }
// buffer = E'F. This computation is done by iterating over the
// f_blocks for each row in the chunk.
@@ -561,6 +565,10 @@
const CompressedRowBlockStructure* bs = A->block_structure();
const double* values = A->values();
for (; row_block_counter < bs->rows.size(); ++row_block_counter) {
+ NoEBlockRowOuterProduct(A, row_block_counter, lhs);
+ if (!rhs) {
+ continue;
+ }
const CompressedRow& row = bs->rows[row_block_counter];
for (int c = 0; c < row.cells.size(); ++c) {
const int block_id = row.cells[c].block_id;
@@ -571,7 +579,6 @@
b + row.block.position,
rhs + lhs_row_layout_[block]);
}
- NoEBlockRowOuterProduct(A, row_block_counter, lhs);
}
}
diff --git a/internal/ceres/schur_jacobi_preconditioner.cc b/internal/ceres/schur_jacobi_preconditioner.cc
index 3ad5dd7..3a3dffd 100644
--- a/internal/ceres/schur_jacobi_preconditioner.cc
+++ b/internal/ceres/schur_jacobi_preconditioner.cc
@@ -88,19 +88,8 @@
const int num_rows = m_->num_rows();
CHECK_GT(num_rows, 0);
- // We need a dummy rhs vector and a dummy b vector since the Schur
- // eliminator combines the computation of the reduced camera matrix
- // with the computation of the right hand side of that linear
- // system.
- //
- // TODO(sameeragarwal): Perhaps its worth refactoring the
- // SchurEliminator::Eliminate function to allow NULL for the rhs. As
- // of now it does not seem to be worth the effort.
- Vector rhs = Vector::Zero(m_->num_rows());
- Vector b = Vector::Zero(A.num_rows());
-
// Compute a subset of the entries of the Schur complement.
- eliminator_->Eliminate(&A, b.data(), D, m_.get(), rhs.data());
+ eliminator_->Eliminate(&A, nullptr, D, m_.get(), nullptr);
m_->Invert();
return true;
}
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index 7cd9bd4..5eebb69 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -338,19 +338,8 @@
const int num_rows = m_->num_rows();
CHECK_GT(num_rows, 0);
- // We need a dummy rhs vector and a dummy b vector since the Schur
- // eliminator combines the computation of the reduced camera matrix
- // with the computation of the right hand side of that linear
- // system.
- //
- // TODO(sameeragarwal): Perhaps its worth refactoring the
- // SchurEliminator::Eliminate function to allow NULL for the rhs. As
- // of now it does not seem to be worth the effort.
- Vector rhs = Vector::Zero(m_->num_rows());
- Vector b = Vector::Zero(A.num_rows());
-
// Compute a subset of the entries of the Schur complement.
- eliminator_->Eliminate(&A, b.data(), D, m_.get(), rhs.data());
+ eliminator_->Eliminate(&A, nullptr, D, m_.get(), nullptr);
// Try factorizing the matrix. For CLUSTER_JACOBI, this should
// always succeed modulo some numerical/conditioning problems. For
