diff --git a/CMakeLists.txt b/CMakeLists.txt index aa444ba31ef6aca9b08fa0400c4201847efd2ac0..d7d567bdb0efdab5f5ae81f7e2961e983d669fca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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/LLT.hpp include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp ) diff --git a/include/eigenpy/decompositions/LLT.hpp b/include/eigenpy/decompositions/LLT.hpp new file mode 100644 index 0000000000000000000000000000000000000000..b3879230767570702d2b807d32a3143b0e30a65b --- /dev/null +++ b/include/eigenpy/decompositions/LLT.hpp @@ -0,0 +1,102 @@ +/* + * Copyright 2020 INRIA + */ + +#ifndef __eigenpy_decomposition_llt_hpp__ +#define __eigenpy_decomposition_llt_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> +#include <Eigen/Cholesky> + +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy +{ + + template<typename _MatrixType> + struct LLTSolverVisitor + : public boost::python::def_visitor< LLTSolverVisitor<_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::LLT<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 LLT factorization from a given matrix.")) + + .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("matrixLLT",&Solver::matrixLLT,bp::arg("self"), + "Returns the LLT 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 LLT 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.") + ; + } + + static void expose() + { + static const std::string classname = "LLT" + scalar_name<Scalar>::shortname(); + expose(classname); + } + + static void expose(const std::string & name) + { + namespace bp = boost::python; + bp::class_<Solver>(name.c_str(), + "Standard Cholesky decomposition (LL^T) of a matrix and associated features.\n\n" + "This class performs a LL^T Cholesky decomposition of a symmetric, positive definite matrix A such that A = LL^* = U^*U, where L is lower triangular.\n\n" + "While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b, for that purpose, we recommend the Cholesky decomposition without square root which is more stable and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other situations like generalised eigen problems with hermitian matrices.", + bp::no_init) + .def(LLTSolverVisitor()); + } + + private: + + static MatrixType matrixL(const Solver & self) { return self.matrixL(); } + static MatrixType matrixU(const Solver & self) { return self.matrixU(); } + + template<typename VectorType> + static VectorType solve(const Solver & self, const VectorType & vec) + { + return self.solve(vec); + } + }; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decomposition_llt_hpp__ diff --git a/src/decompositions/decompositions.cpp b/src/decompositions/decompositions.cpp index 8c5862ffd4c0591064b8fa70f85991c81dc1a5d7..d6d5708691b4897ddd54e4c21d44f78e325e6dfd 100644 --- a/src/decompositions/decompositions.cpp +++ b/src/decompositions/decompositions.cpp @@ -7,6 +7,7 @@ #include "eigenpy/decompositions/EigenSolver.hpp" #include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp" +#include "eigenpy/decompositions/LLT.hpp" namespace eigenpy { @@ -17,6 +18,7 @@ namespace eigenpy EigenSolverVisitor<Eigen::MatrixXd>::expose("EigenSolver"); SelfAdjointEigenSolverVisitor<Eigen::MatrixXd>::expose("SelfAdjointEigenSolver"); + LLTSolverVisitor<Eigen::MatrixXd>::expose("LLT"); { using namespace Eigen;