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/include/ceres/types.h b/include/ceres/types.h index 2c0f0e5..5512340 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h
@@ -101,8 +101,7 @@ JACOBI, // Block diagonal of the Schur complement. This preconditioner may - // only be used with the ITERATIVE_SCHUR solver. Requires - // SuiteSparse/CHOLMOD. + // only be used with the ITERATIVE_SCHUR solver. SCHUR_JACOBI, // Visibility clustering based preconditioners.
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;
diff --git a/jni/Android.mk b/jni/Android.mk index 039dd75..fddb81f 100644 --- a/jni/Android.mk +++ b/jni/Android.mk
@@ -137,6 +137,7 @@ $(CERES_SRC_PATH)/parameter_block_ordering.cc \ $(CERES_SRC_PATH)/partitioned_matrix_view.cc \ $(CERES_SRC_PATH)/polynomial.cc \ + $(CERES_SRC_PATH)/preconditioner.cc \ $(CERES_SRC_PATH)/problem.cc \ $(CERES_SRC_PATH)/problem_impl.cc \ $(CERES_SRC_PATH)/program.cc \ @@ -145,6 +146,7 @@ $(CERES_SRC_PATH)/runtime_numeric_diff_cost_function.cc \ $(CERES_SRC_PATH)/schur_complement_solver.cc \ $(CERES_SRC_PATH)/schur_eliminator.cc \ + $(CERES_SRC_PATH)/schur_jacobi_preconditioner.cc \ $(CERES_SRC_PATH)/scratch_evaluate_preparer.cc \ $(CERES_SRC_PATH)/solver.cc \ $(CERES_SRC_PATH)/solver_impl.cc \