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