Squashed 'third_party/ceres/' content from commit e51e9b4
Change-Id: I763587619d57e594d3fa158dc3a7fe0b89a1743b
git-subtree-dir: third_party/ceres
git-subtree-split: e51e9b46f6ca88ab8b2266d0e362771db6d98067
diff --git a/internal/ceres/subset_preconditioner_test.cc b/internal/ceres/subset_preconditioner_test.cc
new file mode 100644
index 0000000..0285680
--- /dev/null
+++ b/internal/ceres/subset_preconditioner_test.cc
@@ -0,0 +1,202 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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 <memory>
+#include "ceres/subset_preconditioner.h"
+#include "Eigen/Dense"
+#include "Eigen/SparseCore"
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/inner_product_computer.h"
+#include "ceres/internal/eigen.h"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+// TODO(sameeragarwal): Refactor the following two functions out of
+// here and sparse_cholesky_test.cc into a more suitable place.
+template <int UpLoType>
+bool SolveLinearSystemUsingEigen(const Matrix& lhs,
+ const Vector rhs,
+ Vector* solution) {
+ Eigen::LLT<Matrix, UpLoType> llt = lhs.selfadjointView<UpLoType>().llt();
+ if (llt.info() != Eigen::Success) {
+ return false;
+ }
+ *solution = llt.solve(rhs);
+ return (llt.info() == Eigen::Success);
+}
+
+// Use Eigen's Dense Cholesky solver to compute the solution to a
+// sparse linear system.
+bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
+ const Vector& rhs,
+ Vector* solution) {
+ Matrix dense_triangular_lhs;
+ lhs.ToDenseMatrix(&dense_triangular_lhs);
+ if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+ Matrix full_lhs = dense_triangular_lhs.selfadjointView<Eigen::Upper>();
+ return SolveLinearSystemUsingEigen<Eigen::Upper>(full_lhs, rhs, solution);
+ }
+ return SolveLinearSystemUsingEigen<Eigen::Lower>(
+ dense_triangular_lhs, rhs, solution);
+}
+
+typedef ::testing::tuple<SparseLinearAlgebraLibraryType, bool> Param;
+
+std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+ Param param = info.param;
+ std::stringstream ss;
+ ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
+ << (::testing::get<1>(param) ? "Diagonal" : "NoDiagonal");
+ return ss.str();
+}
+
+class SubsetPreconditionerTest : public ::testing::TestWithParam<Param> {
+ protected:
+ virtual void SetUp() {
+ BlockSparseMatrix::RandomMatrixOptions options;
+ options.num_col_blocks = 4;
+ options.min_col_block_size = 1;
+ options.max_col_block_size = 4;
+ options.num_row_blocks = 8;
+ options.min_row_block_size = 1;
+ options.max_row_block_size = 4;
+ options.block_density = 0.9;
+
+ m_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
+ start_row_block_ = m_->block_structure()->rows.size();
+
+ // Ensure that the bottom part of the matrix has the same column
+ // block structure.
+ options.col_blocks = m_->block_structure()->cols;
+ b_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
+ m_->AppendRows(*b_);
+
+ // Create a Identity block diagonal matrix with the same column
+ // block structure.
+ diagonal_ = Vector::Ones(m_->num_cols());
+ block_diagonal_.reset(BlockSparseMatrix::CreateDiagonalMatrix(
+ diagonal_.data(), b_->block_structure()->cols));
+
+ // Unconditionally add the block diagonal to the matrix b_,
+ // because either it is either part of b_ to make it full rank, or
+ // we pass the same diagonal matrix later as the parameter D. In
+ // either case the preconditioner matrix is b_' b + D'D.
+ b_->AppendRows(*block_diagonal_);
+ inner_product_computer_.reset(InnerProductComputer::Create(
+ *b_, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+ inner_product_computer_->Compute();
+ }
+
+ std::unique_ptr<BlockSparseMatrix> m_;
+ std::unique_ptr<BlockSparseMatrix> b_;
+ std::unique_ptr<BlockSparseMatrix> block_diagonal_;
+ std::unique_ptr<InnerProductComputer> inner_product_computer_;
+ std::unique_ptr<Preconditioner> preconditioner_;
+ Vector diagonal_;
+ int start_row_block_;
+};
+
+TEST_P(SubsetPreconditionerTest, foo) {
+ Param param = GetParam();
+ Preconditioner::Options options;
+ options.subset_preconditioner_start_row_block = start_row_block_;
+ options.sparse_linear_algebra_library_type = ::testing::get<0>(param);
+ preconditioner_.reset(new SubsetPreconditioner(options, *m_));
+
+ const bool with_diagonal = ::testing::get<1>(param);
+ if (!with_diagonal) {
+ m_->AppendRows(*block_diagonal_);
+ }
+
+ EXPECT_TRUE(
+ preconditioner_->Update(*m_, with_diagonal ? diagonal_.data() : NULL));
+
+ // Repeatedly apply the preconditioner to random vectors and check
+ // that the preconditioned value is the same as one obtained by
+ // solving the linear system directly.
+ for (int i = 0; i < 5; ++i) {
+ CompressedRowSparseMatrix* lhs = inner_product_computer_->mutable_result();
+ Vector rhs = Vector::Random(lhs->num_rows());
+ Vector expected(lhs->num_rows());
+ EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));
+
+ Vector actual(lhs->num_rows());
+ preconditioner_->RightMultiply(rhs.data(), actual.data());
+
+ Matrix eigen_lhs;
+ lhs->ToDenseMatrix(&eigen_lhs);
+ EXPECT_NEAR((actual - expected).norm() / actual.norm(),
+ 0.0,
+ std::numeric_limits<double>::epsilon() * 10)
+ << "\n"
+ << eigen_lhs << "\n"
+ << expected.transpose() << "\n"
+ << actual.transpose();
+ }
+}
+
+#ifndef CERES_NO_SUITESPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithSuiteSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(SUITE_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_CXSPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithCXSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(CX_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_ACCELERATE_SPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithAccelerateSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(ACCELERATE_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+#ifdef CERES_USE_EIGEN_SPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithEigenSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(EIGEN_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+} // namespace internal
+} // namespace ceres