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
path: root/extern/ceres/internal/ceres/schur_complement_solver.cc
diff options
context:
space:
mode:
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;
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-26 06:01:45 +0300