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

internal/ceres/schur_eliminator.h

@@ -218,9 +218,6 @@ class SchurEliminatorBase {
 // 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 >

internal/ceres/schur_eliminator_impl.h

@@ -176,7 +176,9 @@ Eliminate(const BlockSparseMatrix* A,
           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 @@ Eliminate(const BlockSparseMatrix* A,
         //
         //   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 @@ ChunkDiagonalBlockAndGradient(
             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 @@ NoEBlockRowsUpdate(const BlockSparseMatrix* A,
   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 @@ NoEBlockRowsUpdate(const BlockSparseMatrix* A,
           b + row.block.position,
           rhs + lhs_row_layout_[block]);
     }
-    NoEBlockRowOuterProduct(A, row_block_counter, lhs);
   }
 }
 

internal/ceres/schur_jacobi_preconditioner.cc

@@ -88,19 +88,8 @@ bool SchurJacobiPreconditioner::UpdateImpl(const BlockSparseMatrix& A,
   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;
 }

internal/ceres/visibility_based_preconditioner.cc

@@ -338,19 +338,8 @@ bool VisibilityBasedPreconditioner::UpdateImpl(const BlockSparseMatrix& A,
   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