decompositions: expose Simplicial{LLT,LDLT}

CMakeLists.txt

@@ -119,6 +119,9 @@ set(${PROJECT_NAME}_SOLVERS_HEADERS
 
 set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     include/eigenpy/decompositions/decompositions.hpp
+    include/eigenpy/decompositions/sparse/LLT.hpp
+    include/eigenpy/decompositions/sparse/LDLT.hpp
+    include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
     include/eigenpy/decompositions/EigenSolver.hpp
     include/eigenpy/decompositions/LDLT.hpp
     include/eigenpy/decompositions/LLT.hpp

include/eigenpy/decompositions/sparse/LDLT.hpp

+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_ldlt_hpp__
+#define __eigenpy_decomposition_sparse_ldlt_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+namespace eigenpy {
+
+template <typename _MatrixType, int _UpLo = Eigen::Lower,
+          typename _Ordering =
+              Eigen::AMDOrdering<typename _MatrixType::StorageIndex> >
+struct SimplicialLDLTVisitor
+    : public boost::python::def_visitor<
+          SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering> > {
+  typedef SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor;
+  typedef _MatrixType MatrixType;
+
+  typedef Eigen::SimplicialLDLT<MatrixType> Solver;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+      DenseVectorXs;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+                        MatrixType::Options>
+      DenseMatrixXs;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
+        .def(bp::init<MatrixType>(bp::args("self", "matrix"),
+                                  "Constructs and performs the LDLT "
+                                  "factorization from a given matrix."))
+
+        .def("vectorD", &vectorD, bp::arg("self"),
+             "Returns the diagonal vector D.")
+        .def(SimplicialCholeskyVisitor<Solver>());
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "SimplicialLDLT_" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver, boost::noncopyable>(
+        name.c_str(),
+        "A direct sparse LDLT Cholesky factorizations.\n\n"
+        "This class provides a LDL^T Cholesky factorizations of sparse "
+        "matrices that are selfadjoint and positive definite."
+        "The factorization allows for solving A.X = B where X and B can be "
+        "either dense or sparse.\n\n"
+        "In order to reduce the fill-in, a symmetric permutation P is applied "
+        "prior to the factorization such that the factorized matrix is P A "
+        "P^-1.",
+        bp::no_init)
+        .def(SimplicialLDLTVisitor());
+  }
+
+ private:
+  static DenseVectorXs vectorD(const Solver &self) { return self.vectorD(); }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_ldlt_hpp__

include/eigenpy/decompositions/sparse/LLT.hpp

+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_llt_hpp__
+#define __eigenpy_decomposition_sparse_llt_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+namespace eigenpy {
+
+template <typename _MatrixType, int _UpLo = Eigen::Lower,
+          typename _Ordering =
+              Eigen::AMDOrdering<typename _MatrixType::StorageIndex> >
+struct SimplicialLLTVisitor
+    : public boost::python::def_visitor<
+          SimplicialLLTVisitor<_MatrixType, _UpLo, _Ordering> > {
+  typedef SimplicialLLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor;
+  typedef _MatrixType MatrixType;
+
+  typedef Eigen::SimplicialLLT<MatrixType> Solver;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+      DenseVectorXs;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+                        MatrixType::Options>
+      DenseMatrixXs;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
+        .def(bp::init<MatrixType>(bp::args("self", "matrix"),
+                                  "Constructs and performs the LLT "
+                                  "factorization from a given matrix."))
+
+        .def(SimplicialCholeskyVisitor<Solver>());
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "SimplicialLLT_" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver, boost::noncopyable>(
+        name.c_str(),
+        "A direct sparse LLT Cholesky factorizations.\n\n"
+        "This class provides a LL^T Cholesky factorizations of sparse matrices "
+        "that are selfadjoint and positive definite."
+        "The factorization allows for solving A.X = B where X and B can be "
+        "either dense or sparse.\n\n"
+        "In order to reduce the fill-in, a symmetric permutation P is applied "
+        "prior to the factorization such that the factorized matrix is P A "
+        "P^-1.",
+        bp::no_init)
+        .def(SimplicialLLTVisitor());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_llt_hpp__

include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp

+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+#define __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+
+#include <Eigen/SparseCholesky>
+
+namespace eigenpy {
+
+template <typename SimplicialDerived>
+struct SimplicialCholeskyVisitor
+    : public boost::python::def_visitor<
+          SimplicialCholeskyVisitor<SimplicialDerived> > {
+  typedef SimplicialDerived Solver;
+
+  typedef typename SimplicialDerived::MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+      DenseVectorXs;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+                        MatrixType::Options>
+      DenseMatrixXs;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def("analyzePattern", &Solver::analyzePattern,
+           bp::args("self", "matrix"),
+           "Performs a symbolic decomposition on the sparcity of matrix.\n"
+           "This function is particularly useful when solving for several "
+           "problems having the same structure.")
+
+        .def("cols", &Solver::cols, bp::arg("self"),
+             "Returns the number of columns.")
+        .def("rows", &Solver::rows, bp::arg("self"),
+             "Returns the number of rows.")
+
+        .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("compute",
+             (Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute,
+             bp::args("self", "matrix"),
+             "Computes the sparse Cholesky decomposition of a given matrix.",
+             bp::return_self<>())
+
+        .def("determinant", &Solver::determinant, bp::arg("self"),
+             "Returns the determinant of the underlying matrix from the "
+             "current factorization.")
+
+        .def("factorize", &Solver::factorize, bp::args("self", "matrix"),
+             "Performs a numeric decomposition of a given matrix.\n"
+             "The given matrix must has the same sparcity than the matrix on "
+             "which the symbolic decomposition has been performed.\n"
+             "See also analyzePattern().")
+
+        .def("info", &Solver::info, bp::arg("self"),
+             "NumericalIssue if the input contains INF or NaN values or "
+             "overflow occured. Returns Success otherwise.")
+
+        .def("setShift", &Solver::setShift,
+             (bp::args("self", "offset"), bp::arg("scale") = RealScalar(1)),
+             "Sets the shift parameters that will be used to adjust the "
+             "diagonal coefficients during the numerical factorization.\n"
+             "During the numerical factorization, the diagonal coefficients "
+             "are transformed by the following linear model: d_ii = offset + "
+             "scale * d_ii.\n"
+             "The default is the identity transformation with offset=0, and "
+             "scale=1.",
+             bp::return_self<>())
+
+        .def("solve", &solve<DenseVectorXs>, bp::args("self", "b"),
+             "Returns the solution x of A x = b using the current "
+             "decomposition of A.")
+        .def("solve", &solve<DenseMatrixXs>, bp::args("self", "B"),
+             "Returns the solution X of A X = B using the current "
+             "decomposition of A where B is a right hand side matrix.")
+
+        .def("solve", &solve<MatrixType>, bp::args("self", "B"),
+             "Returns the solution X of A X = B using the current "
+             "decomposition of A where B is a right hand side matrix.");
+  }
+
+ private:
+  static MatrixType matrixL(const Solver &self) { return self.matrixL(); }
+  static MatrixType matrixU(const Solver &self) { return self.matrixU(); }
+
+  template <typename MatrixOrVector>
+  static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+    return self.solve(vec);
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__

src/decompositions/decompositions.cpp

 /*
- * Copyright 2020-2021 INRIA
+ * Copyright 2020-2024 INRIA
  */
 
 #include "eigenpy/decompositions/decompositions.hpp"
@@ -7,6 +7,8 @@
 #include "eigenpy/decompositions/EigenSolver.hpp"
 #include "eigenpy/decompositions/LDLT.hpp"
 #include "eigenpy/decompositions/LLT.hpp"
+#include "eigenpy/decompositions/sparse/LLT.hpp"
+#include "eigenpy/decompositions/sparse/LDLT.hpp"
 #include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp"
 #include "eigenpy/decompositions/minres.hpp"
 #include "eigenpy/fwd.hpp"
@@ -15,6 +17,9 @@ namespace eigenpy {
 void exposeDecompositions() {
   using namespace Eigen;
 
+  typedef Eigen::SparseMatrix<double, Eigen::ColMajor> ColMajorSparseMatrix;
+  //  typedef Eigen::SparseMatrix<double,Eigen::RowMajor> RowMajorSparseMatrix;
+
   EigenSolverVisitor<MatrixXd>::expose("EigenSolver");
   SelfAdjointEigenSolverVisitor<MatrixXd>::expose("SelfAdjointEigenSolver");
   LLTSolverVisitor<MatrixXd>::expose("LLT");
@@ -34,5 +39,11 @@ void exposeDecompositions() {
         .value("ABx_lx", ABx_lx)
         .value("BAx_lx", BAx_lx);
   }
+
+  // Expose sparse decompositions
+  {
+    SimplicialLLTVisitor<ColMajorSparseMatrix>::expose("SimplicialLLT");
+    SimplicialLDLTVisitor<ColMajorSparseMatrix>::expose("SimplicialLDLT");
+  }
 }
 }  // namespace eigenpy