Preconditioner refactoring.
1. Added a Preconditioner interface.
2. SCHUR_JACOBI is now its own class and is independent of
SuiteSparse.
Change-Id: Id912ab19cf3736e61d1b90ddaf5bfba33e877ec4
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index b30f0cc..43e251a 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -72,6 +72,7 @@
parameter_block_ordering.cc
partitioned_matrix_view.cc
polynomial.cc
+ preconditioner.cc
problem.cc
problem_impl.cc
program.cc
@@ -80,6 +81,7 @@
runtime_numeric_diff_cost_function.cc
schur_complement_solver.cc
schur_eliminator.cc
+ schur_jacobi_preconditioner.cc
scratch_evaluate_preparer.cc
solver.cc
solver_impl.cc
diff --git a/internal/ceres/block_jacobi_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc
index 474c37f..5785e5d 100644
--- a/internal/ceres/block_jacobi_preconditioner.cc
+++ b/internal/ceres/block_jacobi_preconditioner.cc
@@ -40,10 +40,9 @@
namespace ceres {
namespace internal {
-BlockJacobiPreconditioner::BlockJacobiPreconditioner(const LinearOperator& A)
+BlockJacobiPreconditioner::BlockJacobiPreconditioner(const BlockSparseMatrixBase& A)
: num_rows_(A.num_rows()),
- block_structure_(
- *(down_cast<const BlockSparseMatrix*>(&A)->block_structure())) {
+ block_structure_(*A.block_structure()) {
// Calculate the amount of storage needed.
int storage_needed = 0;
for (int c = 0; c < block_structure_.cols.size(); ++c) {
@@ -64,11 +63,9 @@
}
}
-BlockJacobiPreconditioner::~BlockJacobiPreconditioner() {
-}
+BlockJacobiPreconditioner::~BlockJacobiPreconditioner() {}
-void BlockJacobiPreconditioner::Update(const LinearOperator& matrix, const double* D) {
- const BlockSparseMatrix& A = *(down_cast<const BlockSparseMatrix*>(&matrix));
+bool BlockJacobiPreconditioner::Update(const BlockSparseMatrixBase& A, const double* D) {
const CompressedRowBlockStructure* bs = A.block_structure();
// Compute the diagonal blocks by block inner products.
@@ -114,6 +111,7 @@
.ldlt()
.solve(Matrix::Identity(size, size));
}
+ return true;
}
void BlockJacobiPreconditioner::RightMultiply(const double* x, double* y) const {
diff --git a/internal/ceres/block_jacobi_preconditioner.h b/internal/ceres/block_jacobi_preconditioner.h
index 51f2655..32221f5 100644
--- a/internal/ceres/block_jacobi_preconditioner.h
+++ b/internal/ceres/block_jacobi_preconditioner.h
@@ -32,38 +32,33 @@
#define CERES_INTERNAL_BLOCK_JACOBI_PRECONDITIONER_H_
#include <vector>
-#include "ceres/linear_operator.h"
+#include "ceres/preconditioner.h"
namespace ceres {
namespace internal {
+class BlockSparseMatrixBase;
struct CompressedRowBlockStructure;
class LinearOperator;
-class SparseMatrix;
-// A block Jacobi preconditioner. This is intended for use with conjugate
-// gradients, or other iterative symmetric solvers. To use the preconditioner,
-// create one by passing a BlockSparseMatrix as the linear operator "A" to the
-// constructor. This fixes the sparsity pattern to the pattern of the matrix
-// A^TA.
+// A block Jacobi preconditioner. This is intended for use with
+// conjugate gradients, or other iterative symmetric solvers. To use
+// the preconditioner, create one by passing a BlockSparseMatrix "A"
+// to the constructor. This fixes the sparsity pattern to the pattern
+// of the matrix A^TA.
//
// Before each use of the preconditioner in a solve with conjugate gradients,
// update the matrix by running Update(A, D). The values of the matrix A are
// inspected to construct the preconditioner. The vector D is applied as the
// D^TD diagonal term.
-class BlockJacobiPreconditioner : public LinearOperator {
+class BlockJacobiPreconditioner : public Preconditioner {
public:
// A must remain valid while the BlockJacobiPreconditioner is.
- BlockJacobiPreconditioner(const LinearOperator& A);
+ BlockJacobiPreconditioner(const BlockSparseMatrixBase& A);
virtual ~BlockJacobiPreconditioner();
- // Update the preconditioner with the values found in A. The sparsity pattern
- // must match that of the A passed to the constructor. D is a vector that
- // must have the same number of rows as A, and is applied as a diagonal in
- // addition to the block diagonals of A.
- void Update(const LinearOperator& A, const double* D);
-
- // LinearOperator interface.
+ // Preconditioner interface
+ virtual bool Update(const BlockSparseMatrixBase& A, const double* D);
virtual void RightMultiply(const double* x, double* y) const;
virtual void LeftMultiply(const double* x, double* y) const;
virtual int num_rows() const { return num_rows_; }
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc
index 5ef521c..e2e799f 100644
--- a/internal/ceres/cgnr_solver.cc
+++ b/internal/ceres/cgnr_solver.cc
@@ -42,11 +42,11 @@
CgnrSolver::CgnrSolver(const LinearSolver::Options& options)
: options_(options),
- jacobi_preconditioner_(NULL) {
+ preconditioner_(NULL) {
}
-LinearSolver::Summary CgnrSolver::Solve(
- LinearOperator* A,
+LinearSolver::Summary CgnrSolver::SolveImpl(
+ BlockSparseMatrixBase* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) {
@@ -60,11 +60,11 @@
// Precondition if necessary.
LinearSolver::PerSolveOptions cg_per_solve_options = per_solve_options;
if (options_.preconditioner_type == JACOBI) {
- if (jacobi_preconditioner_.get() == NULL) {
- jacobi_preconditioner_.reset(new BlockJacobiPreconditioner(*A));
+ if (preconditioner_.get() == NULL) {
+ preconditioner_.reset(new BlockJacobiPreconditioner(*A));
}
- jacobi_preconditioner_->Update(*A, per_solve_options.D);
- cg_per_solve_options.preconditioner = jacobi_preconditioner_.get();
+ preconditioner_->Update(*A, per_solve_options.D);
+ cg_per_solve_options.preconditioner = preconditioner_.get();
} else if (options_.preconditioner_type != IDENTITY) {
LOG(FATAL) << "CGNR only supports IDENTITY and JACOBI preconditioners.";
}
diff --git a/internal/ceres/cgnr_solver.h b/internal/ceres/cgnr_solver.h
index 877b4c4..f61c8d1 100644
--- a/internal/ceres/cgnr_solver.h
+++ b/internal/ceres/cgnr_solver.h
@@ -37,6 +37,8 @@
namespace ceres {
namespace internal {
+class Preconditioner;
+
class BlockJacobiPreconditioner;
// A conjugate gradients on the normal equations solver. This directly solves
@@ -46,17 +48,17 @@
//
// as required for solving for x in the least squares sense. Currently only
// block diagonal preconditioning is supported.
-class CgnrSolver : public LinearSolver {
+class CgnrSolver : public BlockSparseMatrixBaseSolver {
public:
explicit CgnrSolver(const LinearSolver::Options& options);
- virtual Summary Solve(LinearOperator* A,
- const double* b,
- const LinearSolver::PerSolveOptions& per_solve_options,
- double* x);
+ virtual Summary SolveImpl(BlockSparseMatrixBase* A,
+ const double* b,
+ const LinearSolver::PerSolveOptions& per_solve_options,
+ double* x);
private:
const LinearSolver::Options options_;
- scoped_ptr<BlockJacobiPreconditioner> jacobi_preconditioner_;
+ scoped_ptr<Preconditioner> preconditioner_;
CERES_DISALLOW_COPY_AND_ASSIGN(CgnrSolver);
};
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc
index cb2d827..9a66c1e 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -42,6 +42,7 @@
#include "ceres/internal/eigen.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/linear_solver.h"
+#include "ceres/schur_jacobi_preconditioner.h"
#include "ceres/triplet_sparse_matrix.h"
#include "ceres/types.h"
#include "ceres/visibility_based_preconditioner.h"
@@ -93,10 +94,19 @@
cg_per_solve_options.r_tolerance = per_solve_options.r_tolerance;
cg_per_solve_options.q_tolerance = per_solve_options.q_tolerance;
- bool is_preconditioner_good = false;
+ Preconditioner::Options preconditioner_options;
+ preconditioner_options.type = options_.preconditioner_type;
+ preconditioner_options.sparse_linear_algebra_library =
+ options_.sparse_linear_algebra_library;
+ preconditioner_options.use_block_amd = options_.use_block_amd;
+ preconditioner_options.num_threads = options_.num_threads;
+ preconditioner_options.row_block_size = options_.row_block_size;
+ preconditioner_options.e_block_size = options_.e_block_size;
+ preconditioner_options.f_block_size = options_.f_block_size;
+ preconditioner_options.elimination_groups = options_.elimination_groups;
+
switch (options_.preconditioner_type) {
case IDENTITY:
- is_preconditioner_good = true;
break;
case JACOBI:
// We need to strip the constness of the block_diagonal_FtF_inverse
@@ -105,33 +115,41 @@
// that the only method ever called on the preconditioner is the
// RightMultiply which is a const method so we don't need to worry
// about the object getting modified.
- cg_per_solve_options.preconditioner =
- const_cast<BlockSparseMatrix*>(
- schur_complement_->block_diagonal_FtF_inverse());
- is_preconditioner_good = true;
+ preconditioner_.reset(
+ new SparseMatrixPreconditionerWrapper(
+ schur_complement_->block_diagonal_FtF_inverse()));
break;
case SCHUR_JACOBI:
+ if (preconditioner_.get() == NULL) {
+ preconditioner_.reset(
+ new SchurJacobiPreconditioner(*A->block_structure(), preconditioner_options));
+ }
+ break;
case CLUSTER_JACOBI:
case CLUSTER_TRIDIAGONAL:
- if (visibility_based_preconditioner_.get() == NULL) {
- visibility_based_preconditioner_.reset(
- new VisibilityBasedPreconditioner(*A->block_structure(), options_));
+ if (preconditioner_.get() == NULL) {
+ preconditioner_.reset(
+ new VisibilityBasedPreconditioner(*A->block_structure(), preconditioner_options));
}
- is_preconditioner_good =
- visibility_based_preconditioner_->Update(*A, per_solve_options.D);
- cg_per_solve_options.preconditioner =
- visibility_based_preconditioner_.get();
break;
default:
LOG(FATAL) << "Unknown Preconditioner Type";
}
+
+ bool preconditioner_update_was_successful = true;
+ if (preconditioner_.get() != NULL) {
+ preconditioner_update_was_successful =
+ preconditioner_->Update(*A, per_solve_options.D);
+ cg_per_solve_options.preconditioner = preconditioner_.get();
+ }
+
event_logger.AddEvent("Setup");
LinearSolver::Summary cg_summary;
cg_summary.num_iterations = 0;
cg_summary.termination_type = FAILURE;
- if (is_preconditioner_good) {
+ if (preconditioner_update_was_successful) {
cg_summary = cg_solver.Solve(schur_complement_.get(),
schur_complement_->rhs().data(),
cg_per_solve_options,
diff --git a/internal/ceres/iterative_schur_complement_solver.h b/internal/ceres/iterative_schur_complement_solver.h
index 862453c..f8abe04 100644
--- a/internal/ceres/iterative_schur_complement_solver.h
+++ b/internal/ceres/iterative_schur_complement_solver.h
@@ -39,13 +39,12 @@
namespace ceres {
namespace internal {
-class BlockSparseMatrix;
class BlockSparseMatrixBase;
class ImplicitSchurComplement;
-class VisibilityBasedPreconditioner;
+class Preconditioner;
// This class implements an iterative solver for the linear least
-// squares problems that have a bi-partitte sparsity structure common
+// squares problems that have a bi-partite sparsity structure common
// to Structure from Motion problems.
//
// The algorithm used by this solver was developed in a series of
@@ -70,9 +69,7 @@
// "Iterative Methods for Sparse Linear Systems".
class IterativeSchurComplementSolver : public BlockSparseMatrixBaseSolver {
public:
- explicit IterativeSchurComplementSolver(
- const LinearSolver::Options& options);
-
+ explicit IterativeSchurComplementSolver(const LinearSolver::Options& options);
virtual ~IterativeSchurComplementSolver();
private:
@@ -84,7 +81,7 @@
LinearSolver::Options options_;
scoped_ptr<internal::ImplicitSchurComplement> schur_complement_;
- scoped_ptr<VisibilityBasedPreconditioner> visibility_based_preconditioner_;
+ scoped_ptr<Preconditioner> preconditioner_;
Vector reduced_linear_system_solution_;
CERES_DISALLOW_COPY_AND_ASSIGN(IterativeSchurComplementSolver);
};
diff --git a/internal/ceres/preconditioner.cc b/internal/ceres/preconditioner.cc
new file mode 100644
index 0000000..05e539f
--- /dev/null
+++ b/internal/ceres/preconditioner.cc
@@ -0,0 +1,63 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2013 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/preconditioner.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+Preconditioner::~Preconditioner() {
+}
+
+SparseMatrixPreconditionerWrapper::SparseMatrixPreconditionerWrapper(
+ const SparseMatrix* matrix)
+ : matrix_(CHECK_NOTNULL(matrix)) {
+}
+
+SparseMatrixPreconditionerWrapper::~SparseMatrixPreconditionerWrapper() {
+}
+
+bool SparseMatrixPreconditionerWrapper::Update(const BlockSparseMatrixBase& A,
+ const double* D) {
+ return true;
+}
+
+void SparseMatrixPreconditionerWrapper::RightMultiply(const double* x,
+ double* y) const {
+ matrix_->RightMultiply(x, y);
+}
+
+int SparseMatrixPreconditionerWrapper::num_rows() const {
+ return matrix_->num_rows();
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h
new file mode 100644
index 0000000..77280c3
--- /dev/null
+++ b/internal/ceres/preconditioner.h
@@ -0,0 +1,148 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2013 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_PRECONDITIONER_H_
+#define CERES_INTERNAL_PRECONDITIONER_H_
+
+#include <vector>
+#include "ceres/linear_operator.h"
+#include "ceres/sparse_matrix.h"
+
+namespace ceres {
+namespace internal {
+
+class BlockSparseMatrixBase;
+class SparseMatrix;
+
+class Preconditioner : public LinearOperator {
+ public:
+ struct Options {
+ Options()
+ : type(JACOBI),
+ sparse_linear_algebra_library(SUITE_SPARSE),
+ use_block_amd(true),
+ num_threads(1),
+ row_block_size(Dynamic),
+ e_block_size(Dynamic),
+ f_block_size(Dynamic) {
+ }
+
+ PreconditionerType type;
+
+ SparseLinearAlgebraLibraryType sparse_linear_algebra_library;
+
+ // See solver.h for explanation of this option.
+ bool use_block_amd;
+
+ // If possible, how many threads the preconditioner can use.
+ int num_threads;
+
+ // Hints about the order in which the parameter blocks should be
+ // eliminated by the linear solver.
+ //
+ // For example if elimination_groups is a vector of size k, then
+ // the linear solver is informed that it should eliminate the
+ // parameter blocks 0 - elimination_groups[0] - 1 first, and then
+ // elimination_groups[0] - elimination_groups[1] and so on. Within
+ // each elimination group, the linear solver is free to choose how
+ // the parameter blocks are ordered. Different linear solvers have
+ // differing requirements on elimination_groups.
+ //
+ // The most common use is for Schur type solvers, where there
+ // should be at least two elimination groups and the first
+ // elimination group must form an independent set in the normal
+ // equations. The first elimination group corresponds to the
+ // num_eliminate_blocks in the Schur type solvers.
+ vector<int> elimination_groups;
+
+ // If the block sizes in a BlockSparseMatrix are fixed, then in
+ // some cases the Schur complement based solvers can detect and
+ // specialize on them.
+ //
+ // It is expected that these parameters are set programmatically
+ // rather than manually.
+ //
+ // Please see schur_complement_solver.h and schur_eliminator.h for
+ // more details.
+ int row_block_size;
+ int e_block_size;
+ int f_block_size;
+ };
+
+ virtual ~Preconditioner();
+
+ // Update the numerical value of the preconditioner for the linear
+ // system:
+ //
+ // | A | x = |b|
+ // |diag(D)| |0|
+ //
+ // for some vector b. It is important that the matrix A have the
+ // same block structure as the one used to construct this object.
+ //
+ // D can be NULL, in which case its interpreted as a diagonal matrix
+ // of size zero.
+ virtual bool Update(const BlockSparseMatrixBase& A, const double* D) = 0;
+
+ // LinearOperator interface. Since the operator is symmetric,
+ // LeftMultiply and num_cols are just calls to RightMultiply and
+ // num_rows respectively. Update() must be called before
+ // RightMultiply can be called.
+ virtual void RightMultiply(const double* x, double* y) const = 0;
+ virtual void LeftMultiply(const double* x, double* y) const {
+ return RightMultiply(x, y);
+ }
+
+ virtual int num_rows() const = 0;
+ virtual int num_cols() const {
+ return num_rows();
+ }
+};
+
+// Wrap a SparseMatrix object as a preconditioner.
+class SparseMatrixPreconditionerWrapper : public Preconditioner {
+ public:
+ // Wrapper does NOT take ownership of the matrix pointer.
+ SparseMatrixPreconditionerWrapper(const SparseMatrix* matrix);
+ virtual ~SparseMatrixPreconditionerWrapper();
+
+ // Preconditioner interface
+ virtual bool Update(const BlockSparseMatrixBase& A, const double* D);
+ virtual void RightMultiply(const double* x, double* y) const;
+ virtual int num_rows() const;
+
+ private:
+ const SparseMatrix* matrix_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_PRECONDITIONER_H_
diff --git a/internal/ceres/schur_jacobi_preconditioner.cc b/internal/ceres/schur_jacobi_preconditioner.cc
new file mode 100644
index 0000000..33a666e
--- /dev/null
+++ b/internal/ceres/schur_jacobi_preconditioner.cc
@@ -0,0 +1,145 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2013 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/schur_jacobi_preconditioner.h"
+
+#include <utility>
+#include <vector>
+#include "Eigen/Dense"
+#include "ceres/block_random_access_sparse_matrix.h"
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/collections_port.h"
+#include "ceres/detect_structure.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/linear_solver.h"
+#include "ceres/schur_eliminator.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+SchurJacobiPreconditioner::SchurJacobiPreconditioner(
+ const CompressedRowBlockStructure& bs,
+ const Preconditioner::Options& options)
+ : options_(options) {
+ CHECK_GT(options_.elimination_groups.size(), 1);
+ CHECK_GT(options_.elimination_groups[0], 0);
+ const int num_blocks = bs.cols.size() - options_.elimination_groups[0];
+ CHECK_GT(num_blocks, 0)
+ << "Jacobian should have atleast 1 f_block for "
+ << "SCHUR_JACOBI preconditioner.";
+
+ block_size_.resize(num_blocks);
+ set<pair<int, int> > block_pairs;
+
+ int num_block_diagonal_entries = 0;
+ for (int i = 0; i < num_blocks; ++i) {
+ block_size_[i] = bs.cols[i + options_.elimination_groups[0]].size;
+ block_pairs.insert(make_pair(i, i));
+ num_block_diagonal_entries += block_size_[i] * block_size_[i];
+ }
+
+ m_.reset(new BlockRandomAccessSparseMatrix(block_size_, block_pairs));
+ InitEliminator(bs);
+}
+
+SchurJacobiPreconditioner::~SchurJacobiPreconditioner() {
+}
+
+// Initialize the SchurEliminator.
+void SchurJacobiPreconditioner::InitEliminator(
+ const CompressedRowBlockStructure& bs) {
+ LinearSolver::Options eliminator_options;
+
+ eliminator_options.elimination_groups = options_.elimination_groups;
+ eliminator_options.num_threads = options_.num_threads;
+
+ DetectStructure(bs, options_.elimination_groups[0],
+ &eliminator_options.row_block_size,
+ &eliminator_options.e_block_size,
+ &eliminator_options.f_block_size);
+
+ eliminator_.reset(SchurEliminatorBase::Create(eliminator_options));
+ eliminator_->Init(options_.elimination_groups[0], &bs);
+}
+
+// Update the values of the preconditioner matrix and factorize it.
+bool SchurJacobiPreconditioner::Update(const BlockSparseMatrixBase& A,
+ const double* D) {
+ 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());
+ return true;
+}
+
+void SchurJacobiPreconditioner::RightMultiply(const double* x,
+ double* y) const {
+ CHECK_NOTNULL(x);
+ CHECK_NOTNULL(y);
+
+ const double* lhs_values =
+ down_cast<BlockRandomAccessSparseMatrix*>(m_.get())->matrix()->values();
+
+ // This loop can be easily multi-threaded with OpenMP if need be.
+ for (int i = 0; i < block_size_.size(); ++i) {
+ const int block_size = block_size_[i];
+ ConstMatrixRef block(lhs_values, block_size, block_size);
+
+ VectorRef(y, block_size) =
+ block
+ .selfadjointView<Eigen::Upper>()
+ .ldlt()
+ .solve(ConstVectorRef(x, block_size));
+
+ x += block_size;
+ y += block_size;
+ lhs_values += block_size * block_size;
+ }
+}
+
+int SchurJacobiPreconditioner::num_rows() const {
+ return m_->num_rows();
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/schur_jacobi_preconditioner.h b/internal/ceres/schur_jacobi_preconditioner.h
new file mode 100644
index 0000000..3addd73
--- /dev/null
+++ b/internal/ceres/schur_jacobi_preconditioner.h
@@ -0,0 +1,110 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2013 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// Detailed descriptions of these preconditions beyond what is
+// documented here can be found in
+//
+// Bundle Adjustment in the Large
+// S. Agarwal, N. Snavely, S. Seitz & R. Szeliski, ECCV 2010
+// http://www.cs.washington.edu/homes/sagarwal/bal.pdf
+
+#ifndef CERES_INTERNAL_SCHUR_JACOBI_PRECONDITIONER_H_
+#define CERES_INTERNAL_SCHUR_JACOBI_PRECONDITIONER_H_
+
+#include <set>
+#include <vector>
+#include <utility>
+#include "ceres/collections_port.h"
+#include "ceres/internal/macros.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/preconditioner.h"
+
+namespace ceres {
+namespace internal {
+
+class BlockRandomAccessSparseMatrix;
+class BlockSparseMatrixBase;
+struct CompressedRowBlockStructure;
+class SchurEliminatorBase;
+
+// This class implements the SCHUR_JACOBI preconditioner for Structure
+// from Motion/Bundle Adjustment problems. Full mathematical details
+// can be found in
+//
+// Bundle Adjustment in the Large
+// S. Agarwal, N. Snavely, S. Seitz & R. Szeliski, ECCV 2010
+// http://www.cs.washington.edu/homes/sagarwal/bal.pdf
+//
+// Example usage:
+//
+// Preconditioner::Options options;
+// options.preconditioner_type = SCHUR_JACOBI;
+// options.elimination_groups.push_back(num_points);
+// options.elimination_groups.push_back(num_cameras);
+// SchurJacobiPreconditioner preconditioner(
+// *A.block_structure(), options);
+// preconditioner.Update(A, NULL);
+// preconditioner.RightMultiply(x, y);
+//
+class SchurJacobiPreconditioner : public Preconditioner {
+ public:
+ // Initialize the symbolic structure of the preconditioner. bs is
+ // the block structure of the linear system to be solved. It is used
+ // to determine the sparsity structure of the preconditioner matrix.
+ //
+ // It has the same structural requirement as other Schur complement
+ // based solvers. Please see schur_eliminator.h for more details.
+ SchurJacobiPreconditioner(const CompressedRowBlockStructure& bs,
+ const Preconditioner::Options& options);
+ virtual ~SchurJacobiPreconditioner();
+
+ // Preconditioner interface.
+ virtual bool Update(const BlockSparseMatrixBase& A, const double* D);
+ virtual void RightMultiply(const double* x, double* y) const;
+ virtual int num_rows() const;
+
+ private:
+ void InitEliminator(const CompressedRowBlockStructure& bs);
+
+ Preconditioner::Options options_;
+
+ // Sizes of the blocks in the schur complement.
+ vector<int> block_size_;
+ scoped_ptr<SchurEliminatorBase> eliminator_;
+
+ // Preconditioner matrix.
+ scoped_ptr<BlockRandomAccessSparseMatrix> m_;
+ CERES_DISALLOW_COPY_AND_ASSIGN(SchurJacobiPreconditioner);
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_SCHUR_JACOBI_PRECONDITIONER_H_
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index e612c26..6d38535 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -185,8 +185,7 @@
if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
linear_solver_type_used == SPARSE_SCHUR ||
(linear_solver_type_used == ITERATIVE_SCHUR &&
- (preconditioner_type == SCHUR_JACOBI ||
- preconditioner_type == CLUSTER_JACOBI ||
+ (preconditioner_type == CLUSTER_JACOBI ||
preconditioner_type == CLUSTER_TRIDIAGONAL))) {
StringAppendF(&report, "\nSparse Linear Algebra Library %15s\n",
SparseLinearAlgebraLibraryTypeToString(
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 61275fc..c6ee86a 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -1201,12 +1201,6 @@
return NULL;
}
- if (options->preconditioner_type == SCHUR_JACOBI) {
- *error = "SCHUR_JACOBI preconditioner not suppored. Please build Ceres "
- "with SuiteSparse support.";
- return NULL;
- }
-
if (options->preconditioner_type == CLUSTER_JACOBI) {
*error = "CLUSTER_JACOBI preconditioner not suppored. Please build Ceres "
"with SuiteSparse support.";
diff --git a/internal/ceres/system_test.cc b/internal/ceres/system_test.cc
index 4548bd0..597c7ef 100644
--- a/internal/ceres/system_test.cc
+++ b/internal/ceres/system_test.cc
@@ -501,17 +501,19 @@
CONFIGURE(CGNR, SUITE_SPARSE, kAutomaticOrdering, JACOBI);
CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, JACOBI);
+ CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, SCHUR_JACOBI);
#ifndef CERES_NO_SUITESPARSE
- CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, SCHUR_JACOBI);
+
CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, CLUSTER_JACOBI);
CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, CLUSTER_TRIDIAGONAL);
#endif // CERES_NO_SUITESPARSE
CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, JACOBI);
+ CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, SCHUR_JACOBI);
#ifndef CERES_NO_SUITESPARSE
- CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, SCHUR_JACOBI);
+
CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, CLUSTER_JACOBI);
CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, CLUSTER_TRIDIAGONAL);
#endif // CERES_NO_SUITESPARSE
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index ae26d91..a75d6f0 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -65,17 +65,16 @@
#ifndef CERES_NO_SUITESPARSE
VisibilityBasedPreconditioner::VisibilityBasedPreconditioner(
const CompressedRowBlockStructure& bs,
- const LinearSolver::Options& options)
+ const Preconditioner::Options& options)
: options_(options),
num_blocks_(0),
num_clusters_(0),
factor_(NULL) {
CHECK_GT(options_.elimination_groups.size(), 1);
CHECK_GT(options_.elimination_groups[0], 0);
- CHECK(options_.preconditioner_type == SCHUR_JACOBI ||
- options_.preconditioner_type == CLUSTER_JACOBI ||
- options_.preconditioner_type == CLUSTER_TRIDIAGONAL)
- << "Unknown preconditioner type: " << options_.preconditioner_type;
+ CHECK(options_.type == CLUSTER_JACOBI ||
+ options_.type == CLUSTER_TRIDIAGONAL)
+ << "Unknown preconditioner type: " << options_.type;
num_blocks_ = bs.cols.size() - options_.elimination_groups[0];
CHECK_GT(num_blocks_, 0)
<< "Jacobian should have atleast 1 f_block for "
@@ -88,10 +87,7 @@
}
const time_t start_time = time(NULL);
- switch (options_.preconditioner_type) {
- case SCHUR_JACOBI:
- ComputeSchurJacobiSparsity(bs);
- break;
+ switch (options_.type) {
case CLUSTER_JACOBI:
ComputeClusterJacobiSparsity(bs);
break;
@@ -131,24 +127,6 @@
}
}
-// Determine the sparsity structure of the SCHUR_JACOBI
-// preconditioner. SCHUR_JACOBI is an extreme case of a visibility
-// based preconditioner where each camera block corresponds to a
-// cluster and there is no interaction between clusters.
-void VisibilityBasedPreconditioner::ComputeSchurJacobiSparsity(
- const CompressedRowBlockStructure& bs) {
- num_clusters_ = num_blocks_;
- cluster_membership_.resize(num_blocks_);
- cluster_pairs_.clear();
-
- // Each camea block is a member of its own cluster and the only
- // cluster pairs are the self edges (i,i).
- for (int i = 0; i < num_clusters_; ++i) {
- cluster_membership_[i] = i;
- cluster_pairs_.insert(make_pair(i, i));
- }
-}
-
// Determine the sparsity structure of the CLUSTER_JACOBI
// preconditioner. It clusters cameras using their scene
// visibility. The clusters form the diagonal blocks of the
@@ -332,7 +310,6 @@
void VisibilityBasedPreconditioner::InitEliminator(
const CompressedRowBlockStructure& bs) {
LinearSolver::Options eliminator_options;
-
eliminator_options.elimination_groups = options_.elimination_groups;
eliminator_options.num_threads = options_.num_threads;
@@ -366,13 +343,13 @@
// Compute a subset of the entries of the Schur complement.
eliminator_->Eliminate(&A, b.data(), D, m_.get(), rhs.data());
- // Try factorizing the matrix. For SCHUR_JACOBI and CLUSTER_JACOBI,
- // this should always succeed modulo some numerical/conditioning
- // problems. For CLUSTER_TRIDIAGONAL, in general the preconditioner
- // matrix as constructed is not positive definite. However, we will
- // go ahead and try factorizing it. If it works, great, otherwise we
- // scale all the cells in the preconditioner corresponding to the
- // edges in the degree-2 forest and that guarantees positive
+ // Try factorizing the matrix. For CLUSTER_JACOBI, this should
+ // always succeed modulo some numerical/conditioning problems. For
+ // CLUSTER_TRIDIAGONAL, in general the preconditioner matrix as
+ // constructed is not positive definite. However, we will go ahead
+ // and try factorizing it. If it works, great, otherwise we scale
+ // all the cells in the preconditioner corresponding to the edges in
+ // the degree-2 forest and that guarantees positive
// definiteness. The proof of this fact can be found in Lemma 1 in
// "Visibility Based Preconditioning for Bundle Adjustment".
//
@@ -382,10 +359,10 @@
// The scaling only affects the tri-diagonal case, since
// ScaleOffDiagonalBlocks only pays attenion to the cells that
- // belong to the edges of the degree-2 forest. In the SCHUR_JACOBI
- // and the CLUSTER_JACOBI cases, the preconditioner is guaranteed to
- // be positive semidefinite.
- if (!status && options_.preconditioner_type == CLUSTER_TRIDIAGONAL) {
+ // belong to the edges of the degree-2 forest. In the CLUSTER_JACOBI
+ // case, the preconditioner is guaranteed to be positive
+ // semidefinite.
+ if (!status && options_.type == CLUSTER_TRIDIAGONAL) {
VLOG(1) << "Unscaled factorization failed. Retrying with off-diagonal "
<< "scaling";
ScaleOffDiagonalCells();
diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h
index 3246fb8..8a09c78 100644
--- a/internal/ceres/visibility_based_preconditioner.h
+++ b/internal/ceres/visibility_based_preconditioner.h
@@ -29,25 +29,19 @@
// Author: sameeragarwal@google.com (Sameer Agarwal)
//
// Preconditioners for linear systems that arise in Structure from
-// Motion problems. VisibilityBasedPreconditioner implements three
-// preconditioners:
+// Motion problems. VisibilityBasedPreconditioner implements:
//
-// SCHUR_JACOBI
// CLUSTER_JACOBI
// CLUSTER_TRIDIAGONAL
//
// Detailed descriptions of these preconditions beyond what is
// documented here can be found in
//
-// Bundle Adjustment in the Large
-// S. Agarwal, N. Snavely, S. Seitz & R. Szeliski, ECCV 2010
-// http://www.cs.washington.edu/homes/sagarwal/bal.pdf
-//
// Visibility Based Preconditioning for Bundle Adjustment
// A. Kushal & S. Agarwal, submitted to CVPR 2012
// http://www.cs.washington.edu/homes/sagarwal/vbp.pdf
//
-// The three preconditioners share enough code that its most efficient
+// The two preconditioners share enough code that its most efficient
// to implement them as part of the same code base.
#ifndef CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_
@@ -58,11 +52,10 @@
#include <utility>
#include "ceres/collections_port.h"
#include "ceres/graph.h"
-#include "ceres/linear_solver.h"
-#include "ceres/linear_operator.h"
-#include "ceres/suitesparse.h"
#include "ceres/internal/macros.h"
#include "ceres/internal/scoped_ptr.h"
+#include "ceres/preconditioner.h"
+#include "ceres/suitesparse.h"
namespace ceres {
namespace internal {
@@ -72,22 +65,13 @@
struct CompressedRowBlockStructure;
class SchurEliminatorBase;
-// This class implements three preconditioners for Structure from
-// Motion/Bundle Adjustment problems. The name
+// This class implements visibility based preconditioners for
+// Structure from Motion/Bundle Adjustment problems. The name
// VisibilityBasedPreconditioner comes from the fact that the sparsity
// structure of the preconditioner matrix is determined by analyzing
// the visibility structure of the scene, i.e. which cameras see which
// points.
//
-// Strictly speaking, SCHUR_JACOBI is not a visibility based
-// preconditioner but it is an extreme case of CLUSTER_JACOBI, where
-// every cluster contains exactly one camera block. Treating it as a
-// special case of CLUSTER_JACOBI makes it easy to implement as part
-// of the same code base with no significant loss of performance.
-//
-// In the following, we will only discuss CLUSTER_JACOBI and
-// CLUSTER_TRIDIAGONAL.
-//
// The key idea of visibility based preconditioning is to identify
// cameras that we expect have strong interactions, and then using the
// entries in the Schur complement matrix corresponding to these
@@ -130,15 +114,15 @@
//
// LinearSolver::Options options;
// options.preconditioner_type = CLUSTER_JACOBI;
-// options.num_eliminate_blocks = num_points;
+// options.elimination_groups.push_back(num_points);
+// options.elimination_groups.push_back(num_cameras);
// VisibilityBasedPreconditioner preconditioner(
// *A.block_structure(), options);
// preconditioner.Update(A, NULL);
// preconditioner.RightMultiply(x, y);
//
-
#ifndef CERES_NO_SUITESPARSE
-class VisibilityBasedPreconditioner : public LinearOperator {
+class VisibilityBasedPreconditioner : public Preconditioner {
public:
// Initialize the symbolic structure of the preconditioner. bs is
// the block structure of the linear system to be solved. It is used
@@ -146,48 +130,17 @@
//
// It has the same structural requirement as other Schur complement
// based solvers. Please see schur_eliminator.h for more details.
- //
- // LinearSolver::Options::num_eliminate_blocks should be set to the
- // number of e_blocks in the block structure.
- //
- // TODO(sameeragarwal): The use of LinearSolver::Options should
- // ultimately be replaced with Preconditioner::Options and some sort
- // of preconditioner factory along the lines of
- // LinearSolver::CreateLinearSolver. I will wait to do this till I
- // create a general purpose block Jacobi preconditioner for general
- // sparse problems along with a CGLS solver.
VisibilityBasedPreconditioner(const CompressedRowBlockStructure& bs,
- const LinearSolver::Options& options);
+ const Preconditioner::Options& options);
virtual ~VisibilityBasedPreconditioner();
- // Update the numerical value of the preconditioner for the linear
- // system:
- //
- // | A | x = |b|
- // |diag(D)| |0|
- //
- // for some vector b. It is important that the matrix A have the
- // same block structure as the one used to construct this object.
- //
- // D can be NULL, in which case its interpreted as a diagonal matrix
- // of size zero.
- bool Update(const BlockSparseMatrixBase& A, const double* D);
-
-
- // LinearOperator interface. Since the operator is symmetric,
- // LeftMultiply and num_cols are just calls to RightMultiply and
- // num_rows respectively. Update() must be called before
- // RightMultiply can be called.
+ // Preconditioner interface
+ virtual bool Update(const BlockSparseMatrixBase& A, const double* D);
virtual void RightMultiply(const double* x, double* y) const;
- virtual void LeftMultiply(const double* x, double* y) const {
- RightMultiply(x, y);
- }
virtual int num_rows() const;
- virtual int num_cols() const { return num_rows(); }
friend class VisibilityBasedPreconditionerTest;
private:
- void ComputeSchurJacobiSparsity(const CompressedRowBlockStructure& bs);
void ComputeClusterJacobiSparsity(const CompressedRowBlockStructure& bs);
void ComputeClusterTridiagonalSparsity(const CompressedRowBlockStructure& bs);
void InitStorage(const CompressedRowBlockStructure& bs);
@@ -207,7 +160,7 @@
bool IsBlockPairInPreconditioner(int block1, int block2) const;
bool IsBlockPairOffDiagonal(int block1, int block2) const;
- LinearSolver::Options options_;
+ Preconditioner::Options options_;
// Number of parameter blocks in the schur complement.
int num_blocks_;
@@ -249,10 +202,10 @@
#else // SuiteSparse
// If SuiteSparse is not compiled in, the preconditioner is not
// available.
-class VisibilityBasedPreconditioner : public LinearOperator {
+class VisibilityBasedPreconditioner : public Preconditioner {
public:
VisibilityBasedPreconditioner(const CompressedRowBlockStructure& bs,
- const LinearSolver::Options& options) {
+ const Preconditioner::Options& options) {
LOG(FATAL) << "Visibility based preconditioning is not available. Please "
"build Ceres with SuiteSparse.";
}
diff --git a/internal/ceres/visibility_based_preconditioner_test.cc b/internal/ceres/visibility_based_preconditioner_test.cc
index 8c5378d..999024a 100644
--- a/internal/ceres/visibility_based_preconditioner_test.cc
+++ b/internal/ceres/visibility_based_preconditioner_test.cc
@@ -99,7 +99,11 @@
Vector rhs(schur_complement_->num_rows());
scoped_ptr<SchurEliminatorBase> eliminator;
- eliminator.reset(SchurEliminatorBase::Create(options_));
+ LinearSolver::Options eliminator_options;
+ eliminator_options.elimination_groups = options_.elimination_groups;
+ eliminator_options.num_threads = options_.num_threads;
+
+ eliminator.reset(SchurEliminatorBase::Create(eliminator_options));
eliminator->Init(num_eliminate_blocks_, bs);
eliminator->Eliminate(A_.get(), b_.get(), D_.get(),
schur_complement_.get(), rhs.data());
@@ -229,35 +233,14 @@
scoped_array<double> b_;
scoped_array<double> D_;
- LinearSolver::Options options_;
+ Preconditioner::Options options_;
scoped_ptr<VisibilityBasedPreconditioner> preconditioner_;
scoped_ptr<BlockRandomAccessDenseMatrix> schur_complement_;
};
#ifndef CERES_NO_PROTOCOL_BUFFERS
-TEST_F(VisibilityBasedPreconditionerTest, SchurJacobiStructure) {
- options_.preconditioner_type = SCHUR_JACOBI;
- preconditioner_.reset(
- new VisibilityBasedPreconditioner(*A_->block_structure(), options_));
- EXPECT_EQ(get_num_blocks(), num_camera_blocks_);
- EXPECT_EQ(get_num_clusters(), num_camera_blocks_);
- for (int i = 0; i < num_camera_blocks_; ++i) {
- for (int j = 0; j < num_camera_blocks_; ++j) {
- const string msg = StringPrintf("Camera pair: %d %d", i, j);
- SCOPED_TRACE(msg);
- if (i == j) {
- EXPECT_TRUE(IsBlockPairInPreconditioner(i, j));
- EXPECT_FALSE(IsBlockPairOffDiagonal(i, j));
- } else {
- EXPECT_FALSE(IsBlockPairInPreconditioner(i, j));
- EXPECT_TRUE(IsBlockPairOffDiagonal(i, j));
- }
- }
- }
-}
-
TEST_F(VisibilityBasedPreconditionerTest, OneClusterClusterJacobi) {
- options_.preconditioner_type = CLUSTER_JACOBI;
+ options_.type = CLUSTER_JACOBI;
preconditioner_.reset(
new VisibilityBasedPreconditioner(*A_->block_structure(), options_));
@@ -304,7 +287,7 @@
TEST_F(VisibilityBasedPreconditionerTest, ClusterJacobi) {
- options_.preconditioner_type = CLUSTER_JACOBI;
+ options_.type = CLUSTER_JACOBI;
preconditioner_.reset(
new VisibilityBasedPreconditioner(*A_->block_structure(), options_));
@@ -330,7 +313,7 @@
TEST_F(VisibilityBasedPreconditionerTest, ClusterTridiagonal) {
- options_.preconditioner_type = CLUSTER_TRIDIAGONAL;
+ options_.type = CLUSTER_TRIDIAGONAL;
preconditioner_.reset(
new VisibilityBasedPreconditioner(*A_->block_structure(), options_));
static const int kNumClusters = 3;
