diff options
author | Sergey Sharybin <sergey@blender.org> | 2022-05-10 17:36:22 +0300 |
---|---|---|
committer | Sergey Sharybin <sergey@blender.org> | 2022-05-10 18:01:20 +0300 |
commit | 3ad2597a4eca5091031c213445c6583e21097d5f (patch) | |
tree | f909af8ad783d1adea67911ddaf1633ad7f570a9 /extern/ceres/internal/ceres/schur_complement_solver.cc | |
parent | b4b85c5ce2752ea9241cbcfa1ddc3f639ad64262 (diff) |
Update Ceres to latest upstream version 2.1.0temp-ceres_update
This release deprecated the Parameterization API and the new Manifolds
API is to be used instead. This is what was done in the Libmv as part
of this change.
Additionally, remove the bundling scripts. Nowadays those are only
leading to a duplicated work to maintain.
Diffstat (limited to 'extern/ceres/internal/ceres/schur_complement_solver.cc')
-rw-r--r-- | extern/ceres/internal/ceres/schur_complement_solver.cc | 93 |
1 files changed, 42 insertions, 51 deletions
diff --git a/extern/ceres/internal/ceres/schur_complement_solver.cc b/extern/ceres/internal/ceres/schur_complement_solver.cc index 65e7854f9e5..bb442b4280b 100644 --- a/extern/ceres/internal/ceres/schur_complement_solver.cc +++ b/extern/ceres/internal/ceres/schur_complement_solver.cc @@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -46,7 +46,6 @@ #include "ceres/conjugate_gradients_solver.h" #include "ceres/detect_structure.h" #include "ceres/internal/eigen.h" -#include "ceres/lapack.h" #include "ceres/linear_solver.h" #include "ceres/sparse_cholesky.h" #include "ceres/triplet_sparse_matrix.h" @@ -63,14 +62,12 @@ using std::vector; namespace { -class BlockRandomAccessSparseMatrixAdapter : public LinearOperator { +class BlockRandomAccessSparseMatrixAdapter final : public LinearOperator { public: explicit BlockRandomAccessSparseMatrixAdapter( const BlockRandomAccessSparseMatrix& m) : m_(m) {} - virtual ~BlockRandomAccessSparseMatrixAdapter() {} - // y = y + Ax; void RightMultiply(const double* x, double* y) const final { m_.SymmetricRightMultiply(x, y); @@ -88,14 +85,12 @@ class BlockRandomAccessSparseMatrixAdapter : public LinearOperator { const BlockRandomAccessSparseMatrix& m_; }; -class BlockRandomAccessDiagonalMatrixAdapter : public LinearOperator { +class BlockRandomAccessDiagonalMatrixAdapter final : public LinearOperator { public: explicit BlockRandomAccessDiagonalMatrixAdapter( const BlockRandomAccessDiagonalMatrix& m) : m_(m) {} - virtual ~BlockRandomAccessDiagonalMatrixAdapter() {} - // y = y + Ax; void RightMultiply(const double* x, double* y) const final { m_.RightMultiply(x, y); @@ -115,6 +110,14 @@ class BlockRandomAccessDiagonalMatrixAdapter : public LinearOperator { } // namespace +SchurComplementSolver::SchurComplementSolver( + const LinearSolver::Options& options) + : options_(options) { + CHECK_GT(options.elimination_groups.size(), 1); + CHECK_GT(options.elimination_groups[0], 0); + CHECK(options.context != nullptr); +} + LinearSolver::Summary SchurComplementSolver::SolveImpl( BlockSparseMatrix* A, const double* b, @@ -123,7 +126,7 @@ LinearSolver::Summary SchurComplementSolver::SolveImpl( EventLogger event_logger("SchurComplementSolver::Solve"); const CompressedRowBlockStructure* bs = A->block_structure(); - if (eliminator_.get() == NULL) { + if (eliminator_.get() == nullptr) { const int num_eliminate_blocks = options_.elimination_groups[0]; const int num_f_blocks = bs->cols.size() - num_eliminate_blocks; @@ -141,9 +144,9 @@ LinearSolver::Summary SchurComplementSolver::SolveImpl( // mechanism that does not cause binary bloat. if (options_.row_block_size == 2 && options_.e_block_size == 3 && options_.f_block_size == 6 && num_f_blocks == 1) { - eliminator_.reset(new SchurEliminatorForOneFBlock<2, 3, 6>); + eliminator_ = std::make_unique<SchurEliminatorForOneFBlock<2, 3, 6>>(); } else { - eliminator_.reset(SchurEliminatorBase::Create(options_)); + eliminator_ = SchurEliminatorBase::Create(options_); } CHECK(eliminator_); @@ -174,6 +177,12 @@ LinearSolver::Summary SchurComplementSolver::SolveImpl( return summary; } +DenseSchurComplementSolver::DenseSchurComplementSolver( + const LinearSolver::Options& options) + : SchurComplementSolver(options), + cholesky_(DenseCholesky::Create(options)) {} + +DenseSchurComplementSolver::~DenseSchurComplementSolver() = default; // Initialize a BlockRandomAccessDenseMatrix to store the Schur // complement. @@ -187,8 +196,8 @@ void DenseSchurComplementSolver::InitStorage( blocks[j] = bs->cols[i].size; } - set_lhs(new BlockRandomAccessDenseMatrix(blocks)); - set_rhs(new double[lhs()->num_rows()]); + set_lhs(std::make_unique<BlockRandomAccessDenseMatrix>(blocks)); + set_rhs(std::make_unique<double[]>(lhs()->num_rows())); } // Solve the system Sx = r, assuming that the matrix S is stored in a @@ -201,8 +210,7 @@ LinearSolver::Summary DenseSchurComplementSolver::SolveReducedLinearSystem( summary.termination_type = LINEAR_SOLVER_SUCCESS; summary.message = "Success."; - const BlockRandomAccessDenseMatrix* m = - down_cast<const BlockRandomAccessDenseMatrix*>(lhs()); + auto* m = down_cast<BlockRandomAccessDenseMatrix*>(mutable_lhs()); const int num_rows = m->num_rows(); // The case where there are no f blocks, and the system is block @@ -212,26 +220,8 @@ LinearSolver::Summary DenseSchurComplementSolver::SolveReducedLinearSystem( } summary.num_iterations = 1; - - if (options().dense_linear_algebra_library_type == EIGEN) { - Eigen::LLT<Matrix, Eigen::Upper> llt = - ConstMatrixRef(m->values(), num_rows, num_rows) - .selfadjointView<Eigen::Upper>() - .llt(); - if (llt.info() != Eigen::Success) { - summary.termination_type = LINEAR_SOLVER_FAILURE; - summary.message = - "Eigen failure. Unable to perform dense Cholesky factorization."; - return summary; - } - - VectorRef(solution, num_rows) = llt.solve(ConstVectorRef(rhs(), num_rows)); - } else { - VectorRef(solution, num_rows) = ConstVectorRef(rhs(), num_rows); - summary.termination_type = LAPACK::SolveInPlaceUsingCholesky( - num_rows, m->values(), solution, &summary.message); - } - + summary.termination_type = cholesky_->FactorAndSolve( + num_rows, m->mutable_values(), rhs(), solution, &summary.message); return summary; } @@ -243,7 +233,7 @@ SparseSchurComplementSolver::SparseSchurComplementSolver( } } -SparseSchurComplementSolver::~SparseSchurComplementSolver() {} +SparseSchurComplementSolver::~SparseSchurComplementSolver() = default; // Determine the non-zero blocks in the Schur Complement matrix, and // initialize a BlockRandomAccessSparseMatrix object. @@ -303,8 +293,8 @@ void SparseSchurComplementSolver::InitStorage( CHECK_GE(row.cells.front().block_id, num_eliminate_blocks); for (int i = 0; i < row.cells.size(); ++i) { int r_block1_id = row.cells[i].block_id - num_eliminate_blocks; - for (int j = 0; j < row.cells.size(); ++j) { - int r_block2_id = row.cells[j].block_id - num_eliminate_blocks; + for (const auto& cell : row.cells) { + int r_block2_id = cell.block_id - num_eliminate_blocks; if (r_block1_id <= r_block2_id) { block_pairs.insert(make_pair(r_block1_id, r_block2_id)); } @@ -312,8 +302,9 @@ void SparseSchurComplementSolver::InitStorage( } } - set_lhs(new BlockRandomAccessSparseMatrix(blocks_, block_pairs)); - set_rhs(new double[lhs()->num_rows()]); + set_lhs( + std::make_unique<BlockRandomAccessSparseMatrix>(blocks_, block_pairs)); + set_rhs(std::make_unique<double[]>(lhs()->num_rows())); } LinearSolver::Summary SparseSchurComplementSolver::SolveReducedLinearSystem( @@ -338,11 +329,10 @@ LinearSolver::Summary SparseSchurComplementSolver::SolveReducedLinearSystem( const CompressedRowSparseMatrix::StorageType storage_type = sparse_cholesky_->StorageType(); if (storage_type == CompressedRowSparseMatrix::UPPER_TRIANGULAR) { - lhs.reset(CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm)); + lhs = CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm); lhs->set_storage_type(CompressedRowSparseMatrix::UPPER_TRIANGULAR); } else { - lhs.reset( - CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm)); + lhs = CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm); lhs->set_storage_type(CompressedRowSparseMatrix::LOWER_TRIANGULAR); } @@ -373,12 +363,12 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingConjugateGradients( // Only SCHUR_JACOBI is supported over here right now. CHECK_EQ(options().preconditioner_type, SCHUR_JACOBI); - if (preconditioner_.get() == NULL) { - preconditioner_.reset(new BlockRandomAccessDiagonalMatrix(blocks_)); + if (preconditioner_.get() == nullptr) { + preconditioner_ = + std::make_unique<BlockRandomAccessDiagonalMatrix>(blocks_); } - BlockRandomAccessSparseMatrix* sc = down_cast<BlockRandomAccessSparseMatrix*>( - const_cast<BlockRandomAccessMatrix*>(lhs())); + auto* sc = down_cast<BlockRandomAccessSparseMatrix*>(mutable_lhs()); // Extract block diagonal from the Schur complement to construct the // schur_jacobi preconditioner. @@ -404,10 +394,11 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingConjugateGradients( VectorRef(solution, num_rows).setZero(); - std::unique_ptr<LinearOperator> lhs_adapter( - new BlockRandomAccessSparseMatrixAdapter(*sc)); - std::unique_ptr<LinearOperator> preconditioner_adapter( - new BlockRandomAccessDiagonalMatrixAdapter(*preconditioner_)); + std::unique_ptr<LinearOperator> lhs_adapter = + std::make_unique<BlockRandomAccessSparseMatrixAdapter>(*sc); + std::unique_ptr<LinearOperator> preconditioner_adapter = + std::make_unique<BlockRandomAccessDiagonalMatrixAdapter>( + *preconditioner_); LinearSolver::Options cg_options; cg_options.min_num_iterations = options().min_num_iterations; |