Verified Commit 775a1b70 authored by Justin Carpentier's avatar Justin Carpentier
Browse files

decompositions: add LDLT decomposition

parent c775c555
......@@ -90,6 +90,7 @@ SET(${PROJECT_NAME}_SOLVERS_HEADERS
SET(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
include/eigenpy/decompositions/decompositions.hpp
include/eigenpy/decompositions/EigenSolver.hpp
include/eigenpy/decompositions/LDLT.hpp
include/eigenpy/decompositions/LLT.hpp
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
)
......
/*
* Copyright 2020 INRIA
*/
#ifndef __eigenpy_decomposition_ldlt_hpp__
#define __eigenpy_decomposition_ldlt_hpp__
#include <boost/python.hpp>
#include <Eigen/Core>
#include <Eigen/Cholesky>
#include "eigenpy/utils/scalar-name.hpp"
namespace eigenpy
{
template<typename _MatrixType>
struct LDLTSolverVisitor
: public boost::python::def_visitor< LDLTSolverVisitor<_MatrixType> >
{
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1,MatrixType::Options> VectorType;
typedef Eigen::LDLT<MatrixType> Solver;
template<class PyClass>
void visit(PyClass& cl) const
{
namespace bp = boost::python;
cl
.def(bp::init<>("Default constructor"))
.def(bp::init<Eigen::DenseIndex>(bp::arg("size"),
"Default constructor with memory preallocation"))
.def(bp::init<MatrixType>(bp::arg("matrix"),
"Constructs a LDLT factorization from a given matrix."))
.def("isNegative",&Solver::isNegative,bp::arg("self"),
"Returns true if the matrix is negative (semidefinite).")
.def("isPositive",&Solver::isPositive,bp::arg("self"),
"Returns true if the matrix is positive (semidefinite).")
.def("matrixL",&matrixL,bp::arg("self"),
"Returns the lower triangular matrix L.")
.def("matrixU",&matrixU,bp::arg("self"),
"Returns the upper triangular matrix U.")
.def("vectorD",&vectorD,bp::arg("self"),
"Returns the coefficients of the diagonal matrix D.")
.def("transpositionsP",&transpositionsP,bp::arg("self"),
"Returns the permutation matrix P.")
.def("matrixLDLT",&Solver::matrixLDLT,bp::arg("self"),
"Returns the LDLT decomposition matrix.",
bp::return_value_policy<bp::return_by_value>())
.def("rankUpdate",(Solver (Solver::*)(const VectorType &, const RealScalar &))&Solver::template rankUpdate<VectorType>,
bp::args("self","vector","sigma"))
.def("adjoint",&Solver::adjoint,bp::arg("self"),
"Returns the adjoint, that is, a reference to the decomposition itself as if the underlying matrix is self-adjoint.",
bp::return_value_policy<bp::reference_existing_object>())
.def("compute",(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::compute,
bp::args("self","matrix"),
"Computes the LDLT of given matrix.",
bp::return_value_policy<bp::reference_existing_object>())
.def("info",&Solver::info,bp::arg("self"),
"NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise.")
.def("rcond",&Solver::rcond,bp::arg("self"),
"Returns an estimate of the reciprocal condition number of the matrix.")
.def("reconstructedMatrix",&Solver::reconstructedMatrix,bp::arg("self"),
"Returns the matrix represented by the decomposition, i.e., it returns the product: L L^*. This function is provided for debug purpose.")
.def("solve",&solve<VectorType>,bp::args("self","b"),
"Returns the solution x of A x = b using the current decomposition of A.")
.def("setZero",&Solver::setZero,bp::arg("self"),
"Clear any existing decomposition.")
;
}
static void expose()
{
static const std::string classname = "LDLT" + scalar_name<Scalar>::shortname();
expose(classname);
}
static void expose(const std::string & name)
{
namespace bp = boost::python;
bp::class_<Solver>(name.c_str(),
"Robust Cholesky decomposition of a matrix with pivoting.\n\n"
"Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite matrix $ A $ such that $ A = P^TLDL^*P $, where P is a permutation matrix, L is lower triangular with a unit diagonal and D is a diagonal matrix.\n\n"
"The decomposition uses pivoting to ensure stability, so that L will have zeros in the bottom right rank(A) - n submatrix. Avoiding the square root on D also stabilizes the computation.",
bp::no_init)
.def(LDLTSolverVisitor());
}
private:
static MatrixType matrixL(const Solver & self) { return self.matrixL(); }
static MatrixType matrixU(const Solver & self) { return self.matrixU(); }
static VectorType vectorD(const Solver & self) { return self.vectorD(); }
static MatrixType transpositionsP(const Solver & self)
{
return self.transpositionsP() * MatrixType::Identity(self.matrixL().rows(),
self.matrixL().rows());
}
template<typename VectorType>
static VectorType solve(const Solver & self, const VectorType & vec)
{
return self.solve(vec);
}
};
} // namespace eigenpy
#endif // ifndef __eigenpy_decomposition_ldlt_hpp__
......@@ -8,6 +8,7 @@
#include "eigenpy/decompositions/EigenSolver.hpp"
#include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp"
#include "eigenpy/decompositions/LLT.hpp"
#include "eigenpy/decompositions/LDLT.hpp"
namespace eigenpy
{
......@@ -19,6 +20,7 @@ namespace eigenpy
EigenSolverVisitor<Eigen::MatrixXd>::expose("EigenSolver");
SelfAdjointEigenSolverVisitor<Eigen::MatrixXd>::expose("SelfAdjointEigenSolver");
LLTSolverVisitor<Eigen::MatrixXd>::expose("LLT");
LDLTSolverVisitor<Eigen::MatrixXd>::expose("LDLT");
{
using namespace Eigen;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment