Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorSergey Sharybin <sergey@blender.org>2022-05-10 17:36:22 +0300
committerSergey Sharybin <sergey@blender.org>2022-05-10 18:01:20 +0300
commit3ad2597a4eca5091031c213445c6583e21097d5f (patch)
treef909af8ad783d1adea67911ddaf1633ad7f570a9 /extern/ceres/internal/ceres/schur_complement_solver.cc
parentb4b85c5ce2752ea9241cbcfa1ddc3f639ad64262 (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.cc93
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;