saberbud/ceres-solver
1
0
Fork 0
mirror of https://github.com/ceres-solver/ceres-solver.git synced 2026-08-29 16:40:38 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add DenseCholesky

Browse Source 
Like SparseCholesky, the DenseCholesky interface abstracts
away the solution of dense linear systems using Cholesky factorization.
This allows the client code to not worry about the type of dense
linear algebra library being used.

DenseNormalCholeskySolver and DenseSchurComplementSolver code
is considerably simpler as a result.

Change-Id: Ie15f09ee376d5f9a64609e6a55ad83e99c76352a
This commit is contained in:
Sameer Agarwal
2022-01-21 18:33:16 -08:00
parent 77c0c4d09c
commit 6d06e9b98f
11 changed files with 445 additions and 173 deletions
Download Patch File Download Diff File Expand all files Collapse all files
internal/ceres/dense_normal_cholesky_solver.cc
+4 -69
View File
@@ -36,7 +36,6 @@
#include "ceres/blas.h"
#include "ceres/dense_sparse_matrix.h"
#include "ceres/internal/eigen.h"
#include "ceres/lapack.h"
#include "ceres/linear_solver.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
@@ -46,25 +45,13 @@ namespace internal {
DenseNormalCholeskySolver::DenseNormalCholeskySolver(
const LinearSolver::Options& options)
: options_(options) {}
: options_(options), cholesky_(DenseCholesky::Create(options_)) {}
LinearSolver::Summary DenseNormalCholeskySolver::SolveImpl(
DenseSparseMatrix* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) {
if (options_.dense_linear_algebra_library_type == EIGEN) {
return SolveUsingEigen(A, b, per_solve_options, x);
} else {
return SolveUsingLAPACK(A, b, per_solve_options, x);
}
}
LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen(
DenseSparseMatrix* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) {
EventLogger event_logger("DenseNormalCholeskySolver::Solve");
const int num_rows = A->num_rows();
@@ -94,64 +81,12 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen(
LinearSolver::Summary summary;
summary.num_iterations = 1;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
Eigen::LLT<Matrix, Eigen::Upper> llt =
lhs.selfadjointView<Eigen::Upper>().llt();
summary.termination_type = cholesky_->FactorAndSolve(
num_cols, lhs.data(), rhs.data(), x, &summary.message);
event_logger.AddEvent("FactorAndSolve");
if (llt.info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen LLT decomposition failed.";
} else {
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message = "Success.";
}
VectorRef(x, num_cols) = llt.solve(rhs);
event_logger.AddEvent("Solve");
return summary;
}
LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingLAPACK(
DenseSparseMatrix* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) {
EventLogger event_logger("DenseNormalCholeskySolver::Solve");
if (per_solve_options.D != NULL) {
// Temporarily append a diagonal block to the A matrix, but undo
// it before returning the matrix to the user.
A->AppendDiagonal(per_solve_options.D);
}
const int num_cols = A->num_cols();
Matrix lhs(num_cols, num_cols);
event_logger.AddEvent("Setup");
// lhs = A'A
//
// Note: This is a bit delicate, it assumes that the stride on this
// matrix is the same as the number of rows.
BLAS::SymmetricRankKUpdate(
A->num_rows(), num_cols, A->values(), true, 1.0, 0.0, lhs.data());
if (per_solve_options.D != NULL) {
// Undo the modifications to the matrix A.
A->RemoveDiagonal();
}
// TODO(sameeragarwal): Replace this with a gemv call for true blasness.
// rhs = A'b
VectorRef(x, num_cols) =
A->matrix().transpose() * ConstVectorRef(b, A->num_rows());
event_logger.AddEvent("Product");
LinearSolver::Summary summary;
summary.num_iterations = 1;
summary.termination_type = LAPACK::SolveInPlaceUsingCholesky(
num_cols, lhs.data(), x, &summary.message);
event_logger.AddEvent("Solve");
return summary;
}
} // namespace internal
} // namespace ceres