Repository: stack-of-tasks/eigenpy (language: python)

Issue:
Expose remaining Eigen::Decomposition solvers



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8d798ab4b..907cd39e9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - Support for `std::variant` types with `VariantConverter` ([#431](https://github.com/stack-of-tasks/eigenpy/pull/431))
 - Support for `std::unique_ptr` as a return types with `StdUniquePtrCallPolicies` and `boost::python::default_call_policies` ([#433](https://github.com/stack-of-tasks/eigenpy/pull/433))
 - Support for `std::unique_ptr` as an internal reference with `ReturnInternalStdUniquePtr` ([#433](https://github.com/stack-of-tasks/eigenpy/pull/433))
+- Support for `Eigen::Simplicial{LLT,LDLT}` and `Eigen::Cholmod{Simplicial,Supernodal}{LLT,LDLT}` Cholesky de compositions ([#438](https://github.com/stack-of-tasks/eigenpy/pull/438))
 
 ### Fixed
 - Fix the issue of missing exposition of Eigen types with __int64 scalar type ([#426](https://github.com/stack-of-tasks/eigenpy/pull/426))
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e48c9105e..be43048c9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -56,6 +56,15 @@ include("${JRL_CMAKE_MODULES}/apple.cmake")
 option(GENERATE_PYTHON_STUBS
        "Generate the Python stubs associated to the Python library" OFF)
 
+option(BUILD_WITH_CHOLMOD_SUPPORT "Build EigenPy with the Cholmod support" OFF)
+
+if(APPLE)
+  option(BUILD_WITH_ACCELERATE_SUPPORT
+         "Build EigenPy with the Accelerate support" OFF)
+else(APPLE)
+  set(BUILD_WITH_ACCELERATE_SUPPORT FALSE)
+endif(APPLE)
+
 string(REPLACE "-pedantic" "" CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS})
 
 # If needed, fix CMake policy for APPLE systems
@@ -99,6 +108,35 @@ export_boost_default_options()
 find_package(Boost REQUIRED)
 search_for_boost_python(REQUIRED)
 
+if(BUILD_WITH_CHOLMOD_SUPPORT)
+  set(CMAKE_MODULE_PATH ${JRL_CMAKE_MODULES}/find-external/CHOLMOD
+                        ${CMAKE_MODULE_PATH})
+  add_project_dependency(CHOLMOD REQUIRED FIND_EXTERNAL "CHOLMOD")
+  message(
+    STATUS
+      "Build with CHOLDOD support (LGPL). See CHOLMOD/Doc/License.txt for further details."
+  )
+  add_definitions(-DEIGENPY_WITH_CHOLMOD_SUPPORT)
+endif(BUILD_WITH_CHOLMOD_SUPPORT)
+
+if(BUILD_WITH_ACCELERATE_SUPPORT)
+  if(NOT ${Eigen3_VERSION} VERSION_GREATER_EQUAL "3.4.90")
+    message(
+      FATAL_ERROR
+        "Your version of Eigen is too low. Should be at least 3.4.90. Current version is ${Eigen3_VERSION}."
+    )
+  endif()
+
+  set(CMAKE_MODULE_PATH ${JRL_CMAKE_MODULES}/find-external/Accelerate
+                        ${CMAKE_MODULE_PATH})
+  find_package(
+    Accelerate REQUIRED # FIND_EXTERNAL "Accelerate" # We don't export yet as
+                        # there might be an issue on AMR64 platforms
+  )
+  message(STATUS "Build with Accelerate support framework.")
+  add_definitions(-DEIGENPY_WITH_ACCELERATE_SUPPORT)
+endif(BUILD_WITH_ACCELERATE_SUPPORT)
+
 # ----------------------------------------------------
 # --- INCLUDE ----------------------------------------
 # ----------------------------------------------------
@@ -117,9 +155,39 @@ set(${PROJECT_NAME}_SOLVERS_HEADERS
     include/eigenpy/solvers/BasicPreconditioners.hpp
     include/eigenpy/solvers/BFGSPreconditioners.hpp)
 
+set(${PROJECT_NAME}_EIGEN_HEADERS include/eigenpy/eigen/EigenBase.hpp)
+
+set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_CHOLMOD_HEADERS
+    include/eigenpy/decompositions/sparse/cholmod/CholmodBase.hpp
+    include/eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp
+    include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLDLT.hpp
+    include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLLT.hpp
+    include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp)
+
+set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_ACCELERATE_HEADERS
+    include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp)
+
+set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
+    include/eigenpy/decompositions/sparse/LLT.hpp
+    include/eigenpy/decompositions/sparse/LDLT.hpp
+    include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
+    include/eigenpy/decompositions/sparse/SparseSolverBase.hpp)
+
+if(BUILD_WITH_CHOLMOD_SUPPORT)
+  list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
+       ${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_CHOLMOD_HEADERS})
+endif(BUILD_WITH_CHOLMOD_SUPPORT)
+
+if(BUILD_WITH_ACCELERATE_SUPPORT)
+  list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
+       ${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_ACCELERATE_HEADERS})
+endif(BUILD_WITH_ACCELERATE_SUPPORT)
+
 set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
+    ${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS}
     include/eigenpy/decompositions/decompositions.hpp
     include/eigenpy/decompositions/EigenSolver.hpp
+    include/eigenpy/decompositions/PermutationMatrix.hpp
     include/eigenpy/decompositions/LDLT.hpp
     include/eigenpy/decompositions/LLT.hpp
     include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
@@ -128,6 +196,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
 set(${PROJECT_NAME}_HEADERS
     ${${PROJECT_NAME}_UTILS_HEADERS}
     ${${PROJECT_NAME}_SOLVERS_HEADERS}
+    ${${PROJECT_NAME}_EIGEN_HEADERS}
     ${${PROJECT_NAME}_DECOMPOSITIONS_HEADERS}
     include/eigenpy/alignment.hpp
     include/eigenpy/computation-info.hpp
@@ -188,6 +257,16 @@ set(${PROJECT_NAME}_SOLVERS_SOURCES src/solvers/preconditioners.cpp
 set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/decompositions.cpp)
 
+if(BUILD_WITH_CHOLMOD_SUPPORT)
+  list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
+       src/decompositions/cholmod.cpp)
+endif(BUILD_WITH_CHOLMOD_SUPPORT)
+
+if(BUILD_WITH_ACCELERATE_SUPPORT)
+  list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
+       src/decompositions/accelerate.cpp)
+endif(BUILD_WITH_ACCELERATE_SUPPORT)
+
 set(${PROJECT_NAME}_SOURCES
     ${${PROJECT_NAME}_SOLVERS_SOURCES}
     ${${PROJECT_NAME}_DECOMPOSITIONS_SOURCES}
@@ -239,6 +318,19 @@ modernize_target_link_libraries(
   ${NUMPY_INCLUDE_DIRS}
   ${PYTHON_INCLUDE_DIR})
 
+# Links against CholMod
+if(BUILD_WITH_CHOLMOD_SUPPORT)
+  modernize_target_link_libraries(${PROJECT_NAME} SCOPE PUBLIC TARGETS
+                                  CHOLMOD::CHOLMOD)
+endif(BUILD_WITH_CHOLMOD_SUPPORT)
+
+# Links against accelerate
+if(BUILD_WITH_ACCELERATE_SUPPORT)
+  # modernize_target_link_libraries(${PROJECT_NAME} SCOPE PUBLIC TARGETS
+  # Accelerate)
+  target_link_libraries(${PROJECT_NAME} PRIVATE "-framework accelerate")
+endif(BUILD_WITH_ACCELERATE_SUPPORT)
+
 if(SUFFIX_SO_VERSION)
   set_target_properties(${PROJECT_NAME} PROPERTIES SOVERSION ${PROJECT_VERSION})
 endif(SUFFIX_SO_VERSION)
diff --git a/cmake b/cmake
index dd1d28f51..b3a58c457 160000
--- a/cmake
+++ b/cmake
@@ -1,1 +1,1 @@
-Subproject commit dd1d28f510754e85bf802189e29d601116df6699
+Subproject commit b3a58c457f5510a7ae7ece501b4bc2aaf2ae5f87
diff --git a/include/eigenpy/decompositions/LDLT.hpp b/include/eigenpy/decompositions/LDLT.hpp
index 347405411..2a635a217 100644
--- a/include/eigenpy/decompositions/LDLT.hpp
+++ b/include/eigenpy/decompositions/LDLT.hpp
@@ -10,6 +10,7 @@
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/utils/scalar-name.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
 
 namespace eigenpy {
 
@@ -36,6 +37,8 @@ struct LDLTSolverVisitor
             bp::args("self", "matrix"),
             "Constructs a LDLT factorization from a given matrix."))
 
+        .def(EigenBaseVisitor<Solver>())
+
         .def("isNegative", &Solver::isNegative, bp::arg("self"),
              "Returns true if the matrix is negative (semidefinite).")
         .def("isPositive", &Solver::isPositive, bp::arg("self"),
diff --git a/include/eigenpy/decompositions/LLT.hpp b/include/eigenpy/decompositions/LLT.hpp
index b0a1e1b57..1996d6ccb 100644
--- a/include/eigenpy/decompositions/LLT.hpp
+++ b/include/eigenpy/decompositions/LLT.hpp
@@ -10,6 +10,7 @@
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/utils/scalar-name.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
 
 namespace eigenpy {
 
@@ -36,6 +37,8 @@ struct LLTSolverVisitor
             bp::args("self", "matrix"),
             "Constructs a LLT factorization from a given matrix."))
 
+        .def(EigenBaseVisitor<Solver>())
+
         .def("matrixL", &matrixL, bp::arg("self"),
              "Returns the lower triangular matrix L.")
         .def("matrixU", &matrixU, bp::arg("self"),
diff --git a/include/eigenpy/decompositions/PermutationMatrix.hpp b/include/eigenpy/decompositions/PermutationMatrix.hpp
new file mode 100644
index 000000000..48d8ce7b8
--- /dev/null
+++ b/include/eigenpy/decompositions/PermutationMatrix.hpp
@@ -0,0 +1,106 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_permutation_matrix_hpp__
+#define __eigenpy_decomposition_permutation_matrix_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+
+namespace eigenpy {
+
+template <int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime,
+          typename StorageIndex_ = int>
+struct PermutationMatrixVisitor
+    : public boost::python::def_visitor<PermutationMatrixVisitor<
+          SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> > {
+  typedef StorageIndex_ StorageIndex;
+  typedef Eigen::PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime,
+                                   StorageIndex>
+      PermutationMatrix;
+  typedef typename PermutationMatrix::DenseMatrixType DenseMatrixType;
+  typedef PermutationMatrix Self;
+  typedef Eigen::Matrix<StorageIndex, SizeAtCompileTime, 1, 0,
+                        MaxSizeAtCompileTime, 1>
+      VectorIndex;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def(bp::init<const Eigen::DenseIndex>(bp::args("self", "size"),
+                                             "Default constructor"))
+        .def(bp::init<VectorIndex>(
+            bp::args("self", "indices"),
+            "The indices array has the meaning that the permutations sends "
+            "each integer i to indices[i].\n"
+            "It is your responsibility to check that the indices array that "
+            "you passes actually describes a permutation, i.e., each value "
+            "between 0 and n-1 occurs exactly once, where n is the array's "
+            "size."))
+
+        .def(
+            "indices",
+            +[](const PermutationMatrix &self) {
+              return VectorIndex(self.indices());
+            },
+            bp::arg("self"), "The stored array representing the permutation.")
+
+        .def("applyTranspositionOnTheLeft",
+             &PermutationMatrix::applyTranspositionOnTheLeft,
+             bp::args("self", "i", "j"),
+             "Multiplies self by the transposition (ij) on the left.",
+             bp::return_self<>())
+        .def("applyTranspositionOnTheRight",
+             &PermutationMatrix::applyTranspositionOnTheRight,
+             bp::args("self", "i", "j"),
+             "Multiplies self by the transposition (ij) on the right.",
+             bp::return_self<>())
+
+        .def("setIdentity",
+             (void(PermutationMatrix::*)()) & PermutationMatrix::setIdentity,
+             bp::arg("self"),
+             "Sets self to be the identity permutation matrix.")
+        .def("setIdentity",
+             (void(PermutationMatrix::*)(Eigen::DenseIndex)) &
+                 PermutationMatrix::setIdentity,
+             bp::args("self", "size"),
+             "Sets self to be the identity permutation matrix of given size.")
+
+        .def("toDenseMatrix", &PermutationMatrix::toDenseMatrix,
+             bp::arg("self"),
+             "Returns a numpy array object initialized from this permutation "
+             "matrix.")
+
+        .def(
+            "transpose",
+            +[](const PermutationMatrix &self) -> PermutationMatrix {
+              return self.transpose();
+            },
+            bp::arg("self"), "Returns the tranpose permutation matrix.")
+        .def(
+            "inverse",
+            +[](const PermutationMatrix &self) -> PermutationMatrix {
+              return self.inverse();
+            },
+            bp::arg("self"), "Returns the inverse permutation matrix.")
+
+        .def("resize", &PermutationMatrix::resize, bp::args("self", "size"),
+             "Resizes to given size.")
+
+        .def(bp::self * bp::self)
+        .def(EigenBaseVisitor<Self>());
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<PermutationMatrix>(name.c_str(),
+                                  "Permutation matrix.\n"
+                                  "This class represents a permutation matrix, "
+                                  "internally stored as a vector of integers.",
+                                  bp::no_init)
+        .def(PermutationMatrixVisitor());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_permutation_matrix_hpp__
diff --git a/include/eigenpy/decompositions/decompositions.hpp b/include/eigenpy/decompositions/decompositions.hpp
index 8e6668d5c..fd55626a4 100644
--- a/include/eigenpy/decompositions/decompositions.hpp
+++ b/include/eigenpy/decompositions/decompositions.hpp
@@ -9,6 +9,15 @@
 
 namespace eigenpy {
 void EIGENPY_DLLAPI exposeDecompositions();
+
+#ifdef EIGENPY_WITH_CHOLMOD_SUPPORT
+void EIGENPY_DLLAPI exposeCholmod();
+#endif
+
+#ifdef EIGENPY_WITH_ACCELERATE_SUPPORT
+void EIGENPY_DLLAPI exposeAccelerate();
+#endif
+
 }  // namespace eigenpy
 
 #endif  // define __eigenpy_decompositions_decompositions_hpp__
diff --git a/include/eigenpy/decompositions/sparse/LDLT.hpp b/include/eigenpy/decompositions/sparse/LDLT.hpp
new file mode 100644
index 000000000..d282442cb
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/LDLT.hpp
@@ -0,0 +1,71 @@
+/*
+ * 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__
diff --git a/include/eigenpy/decompositions/sparse/LLT.hpp b/include/eigenpy/decompositions/sparse/LLT.hpp
new file mode 100644
index 000000000..199c25734
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/LLT.hpp
@@ -0,0 +1,66 @@
+/*
+ * 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__
diff --git a/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp b/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
new file mode 100644
index 000000000..ad109140b
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
@@ -0,0 +1,94 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+#define __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+#include "eigenpy/decompositions/sparse/SparseSolverBase.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(EigenBaseVisitor<Solver>())
+        .def(SparseSolverBaseVisitor<Solver>())
+
+        .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("permutationP", &Solver::permutationP, bp::arg("self"),
+             "Returns the permutation P.",
+             bp::return_value_policy<bp::copy_const_reference>())
+        .def("permutationPinv", &Solver::permutationPinv, bp::arg("self"),
+             "Returns the inverse P^-1 of the permutation P.",
+             bp::return_value_policy<bp::copy_const_reference>());
+  }
+
+ private:
+  static MatrixType matrixL(const Solver &self) { return self.matrixL(); }
+  static MatrixType matrixU(const Solver &self) { return self.matrixU(); }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
diff --git a/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp b/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
new file mode 100644
index 000000000..2afc1b917
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
@@ -0,0 +1,54 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
+#define __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+
+#include <Eigen/SparseCholesky>
+
+namespace eigenpy {
+
+template <typename SimplicialDerived>
+struct SparseSolverBaseVisitor
+    : public boost::python::def_visitor<
+          SparseSolverBaseVisitor<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("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:
+  template <typename MatrixOrVector>
+  static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+    return self.solve(vec);
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
diff --git a/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp b/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp
new file mode 100644
index 000000000..903c086fb
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp
@@ -0,0 +1,73 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__
+#define __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+#include "eigenpy/decompositions/sparse/SparseSolverBase.hpp"
+
+#include <Eigen/AccelerateSupport>
+
+namespace eigenpy {
+
+template <typename AccelerateDerived>
+struct AccelerateImplVisitor : public boost::python::def_visitor<
+                                   AccelerateImplVisitor<AccelerateDerived> > {
+  typedef AccelerateDerived Solver;
+
+  typedef typename AccelerateDerived::MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef MatrixType CholMatrixType;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+
+  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(EigenBaseVisitor<Solver>())
+        .def(SparseSolverBaseVisitor<Solver>())
+
+        .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("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("setOrder", &Solver::setOrder, bp::arg("self"), "Set order");
+  }
+
+  static void expose(const std::string &name, const std::string &doc = "") {
+    bp::class_<Solver, boost::noncopyable>(name.c_str(), doc.c_str(),
+                                           bp::no_init)
+        .def(AccelerateImplVisitor())
+
+        .def(bp::init<>(bp::arg("self"), "Default constructor"))
+        .def(bp::init<MatrixType>(bp::args("self", "matrix"),
+                                  "Constructs and performs the "
+                                  "factorization from a given matrix."));
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__
diff --git a/include/eigenpy/decompositions/sparse/cholmod/CholmodBase.hpp b/include/eigenpy/decompositions/sparse/cholmod/CholmodBase.hpp
new file mode 100644
index 000000000..9df15099e
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/cholmod/CholmodBase.hpp
@@ -0,0 +1,75 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_cholmod_cholmod_base_hpp__
+#define __eigenpy_decomposition_sparse_cholmod_cholmod_base_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+#include "eigenpy/decompositions/sparse/SparseSolverBase.hpp"
+
+#include <Eigen/CholmodSupport>
+
+namespace eigenpy {
+
+template <typename CholdmodDerived>
+struct CholmodBaseVisitor
+    : public boost::python::def_visitor<CholmodBaseVisitor<CholdmodDerived> > {
+  typedef CholdmodDerived Solver;
+
+  typedef typename CholdmodDerived::MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef MatrixType CholMatrixType;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+
+  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(EigenBaseVisitor<Solver>())
+        .def(SparseSolverBaseVisitor<Solver>())
+
+        .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("logDeterminant", &Solver::logDeterminant, bp::arg("self"),
+             "Returns the log determinant of the underlying matrix from the "
+             "current factorization.")
+
+        .def("setShift", &Solver::setShift, (bp::args("self", "offset")),
+             "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 + "
+             "d_ii.\n"
+             "The default is the identity transformation with offset=0.",
+             bp::return_self<>());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decomposition_sparse_cholmod_cholmod_base_hpp__
diff --git a/include/eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp b/include/eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp
new file mode 100644
index 000000000..0e2076ca9
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp
@@ -0,0 +1,32 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_cholmod_cholmod_decomposition_hpp__
+#define __eigenpy_decomposition_sparse_cholmod_cholmod_decomposition_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/decompositions/sparse/cholmod/CholmodBase.hpp"
+
+namespace eigenpy {
+
+template <typename CholdmodDerived>
+struct CholmodDecompositionVisitor
+    : public boost::python::def_visitor<
+          CholmodDecompositionVisitor<CholdmodDerived> > {
+  typedef CholdmodDerived Solver;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl
+
+        .def(CholmodBaseVisitor<Solver>())
+        .def("setMode", &Solver::setMode, bp::args("self", "mode"),
+             "Set the mode for the Cholesky decomposition.");
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef
+        // __eigenpy_decomposition_sparse_cholmod_cholmod_decomposition_hpp__
diff --git a/include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLDLT.hpp b/include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLDLT.hpp
new file mode 100644
index 000000000..0e257022f
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLDLT.hpp
@@ -0,0 +1,63 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_cholmod_cholmod_simplicial_ldlt_hpp__
+#define __eigenpy_decomposition_sparse_cholmod_cholmod_simplicial_ldlt_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+namespace eigenpy {
+
+template <typename MatrixType_, int UpLo_ = Eigen::Lower>
+struct CholmodSimplicialLDLTVisitor
+    : public boost::python::def_visitor<
+          CholmodSimplicialLDLTVisitor<MatrixType_, UpLo_> > {
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+
+  typedef Eigen::CholmodSimplicialLDLT<MatrixType_, UpLo_> Solver;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl
+
+        .def(CholmodBaseVisitor<Solver>())
+        .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."))
+
+        ;
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "CholmodSimplicialLDLT_" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver, boost::noncopyable>(
+        name.c_str(),
+        "A simplicial direct Cholesky (LDLT) factorization and solver based on "
+        "Cholmod.\n\n"
+        "This class allows to solve for A.X = B sparse linear problems via a "
+        "simplicial LL^T Cholesky factorization using the Cholmod library."
+        "This simplicial variant is equivalent to Eigen's built-in "
+        "SimplicialLDLT class."
+        "Therefore, it has little practical interest. The sparse matrix A must "
+        "be selfadjoint and positive definite."
+        "The vectors or matrices X and B can be either dense or sparse.",
+        bp::no_init)
+        .def(CholmodSimplicialLDLTVisitor());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef
+        // __eigenpy_decomposition_sparse_cholmod_cholmod_simplicial_ldlt_hpp__
diff --git a/include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLLT.hpp b/include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLLT.hpp
new file mode 100644
index 000000000..c84b2a1e6
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLLT.hpp
@@ -0,0 +1,63 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_cholmod_cholmod_simplicial_llt_hpp__
+#define __eigenpy_decomposition_sparse_cholmod_cholmod_simplicial_llt_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+namespace eigenpy {
+
+template <typename MatrixType_, int UpLo_ = Eigen::Lower>
+struct CholmodSimplicialLLTVisitor
+    : public boost::python::def_visitor<
+          CholmodSimplicialLLTVisitor<MatrixType_, UpLo_> > {
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+
+  typedef Eigen::CholmodSimplicialLLT<MatrixType_, UpLo_> Solver;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl
+
+        .def(CholmodBaseVisitor<Solver>())
+        .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."))
+
+        ;
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "CholmodSimplicialLLT_" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver, boost::noncopyable>(
+        name.c_str(),
+        "A simplicial direct Cholesky (LLT) factorization and solver based on "
+        "Cholmod.\n\n"
+        "This class allows to solve for A.X = B sparse linear problems via a "
+        "simplicial LL^T Cholesky factorization using the Cholmod library."
+        "This simplicial variant is equivalent to Eigen's built-in "
+        "SimplicialLLT class."
+        "Therefore, it has little practical interest. The sparse matrix A must "
+        "be selfadjoint and positive definite."
+        "The vectors or matrices X and B can be either dense or sparse.",
+        bp::no_init)
+        .def(CholmodSimplicialLLTVisitor());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef
+        // __eigenpy_decomposition_sparse_cholmod_cholmod_simplicial_llt_hpp__
diff --git a/include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp b/include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp
new file mode 100644
index 000000000..ca8fdfe67
--- /dev/null
+++ b/include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp
@@ -0,0 +1,62 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_sparse_cholmod_cholmod_supernodal_llt_hpp__
+#define __eigenpy_decomposition_sparse_cholmod_cholmod_supernodal_llt_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+namespace eigenpy {
+
+template <typename MatrixType_, int UpLo_ = Eigen::Lower>
+struct CholmodSupernodalLLTVisitor
+    : public boost::python::def_visitor<
+          CholmodSupernodalLLTVisitor<MatrixType_, UpLo_> > {
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+
+  typedef Eigen::CholmodSupernodalLLT<MatrixType_, UpLo_> Solver;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl
+
+        .def(CholmodBaseVisitor<Solver>())
+        .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."))
+
+        ;
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "CholmodSupernodalLLT_" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver, boost::noncopyable>(
+        name.c_str(),
+        "A supernodal direct Cholesky (LLT) factorization and solver based on "
+        "Cholmod.\n\n"
+        "This class allows to solve for A.X = B sparse linear problems via a "
+        "supernodal LL^T Cholesky factorization using the Cholmod library."
+        "This supernodal variant performs best on dense enough problems, e.g., "
+        "3D FEM, or very high order 2D FEM."
+        "The sparse matrix A must be selfadjoint and positive definite. The "
+        "vectors or matrices X and B can be either dense or sparse.",
+        bp::no_init)
+        .def(CholmodSupernodalLLTVisitor());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef
+        // __eigenpy_decomposition_sparse_cholmod_cholmod_supernodal_llt_hpp__
diff --git a/include/eigenpy/eigen-typedef.hpp b/include/eigenpy/eigen-typedef.hpp
index d195d11a6..9d343ff2c 100644
--- a/include/eigenpy/eigen-typedef.hpp
+++ b/include/eigenpy/eigen-typedef.hpp
@@ -32,6 +32,7 @@
   EIGENPY_MAKE_FIXED_TYPEDEFS(Type, Options, TypeSuffix, 2)           \
   EIGENPY_MAKE_FIXED_TYPEDEFS(Type, Options, TypeSuffix, 3)           \
   EIGENPY_MAKE_FIXED_TYPEDEFS(Type, Options, TypeSuffix, 4)           \
-  EIGENPY_MAKE_TYPEDEFS(Type, Options, TypeSuffix, 1, 1)
+  EIGENPY_MAKE_TYPEDEFS(Type, Options, TypeSuffix, 1, 1)              \
+  typedef Eigen::SparseMatrix<Scalar, Options> SparseMatrixX##TypeSuffix
 
 #endif  // ifndef __eigenpy_eigen_typedef_hpp__
diff --git a/include/eigenpy/eigen/EigenBase.hpp b/include/eigenpy/eigen/EigenBase.hpp
new file mode 100644
index 000000000..262e7b2d9
--- /dev/null
+++ b/include/eigenpy/eigen/EigenBase.hpp
@@ -0,0 +1,28 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_eigen_eigen_base_hpp__
+#define __eigenpy_eigen_eigen_base_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+
+namespace eigenpy {
+
+template <typename Derived>
+struct EigenBaseVisitor
+    : public boost::python::def_visitor<EigenBaseVisitor<Derived> > {
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def("cols", &Derived::cols, bp::arg("self"),
+           "Returns the number of columns.")
+        .def("rows", &Derived::rows, bp::arg("self"),
+             "Returns the number of rows.")
+        .def("size", &Derived::rows, bp::arg("self"),
+             "Returns the number of coefficients, which is rows()*cols().");
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_eigen_eigen_base_hpp__
diff --git a/include/eigenpy/eigenpy.hpp b/include/eigenpy/eigenpy.hpp
index 97407d349..5bbdc1b63 100644
--- a/include/eigenpy/eigenpy.hpp
+++ b/include/eigenpy/eigenpy.hpp
@@ -58,7 +58,7 @@ EIGEN_DONT_INLINE void exposeType() {
   ENABLE_SPECIFIC_MATRIX_TYPE(RowVectorXs);
   ENABLE_SPECIFIC_MATRIX_TYPE(MatrixXs);
 
-  enableEigenPySpecific<Eigen::SparseMatrix<Scalar, Options> >();
+  enableEigenPySpecific<SparseMatrixXs>();
 }
 
 template <typename Scalar>
diff --git a/include/eigenpy/utils/is-approx.hpp b/include/eigenpy/utils/is-approx.hpp
index 989c98bc3..45c178b19 100644
--- a/include/eigenpy/utils/is-approx.hpp
+++ b/include/eigenpy/utils/is-approx.hpp
@@ -1,29 +1,45 @@
 /*
- * Copyright 2020 INRIA
+ * Copyright 2020-2024 INRIA
  */
 
-#ifndef __eigenpy_utils_scalar_is_approx_hpp__
-#define __eigenpy_utils_scalar_is_approx_hpp__
+#ifndef __eigenpy_utils_is_approx_hpp__
+#define __eigenpy_utils_is_approx_hpp__
 
 #include <Eigen/Core>
+#include <Eigen/SparseCore>
 
 namespace eigenpy {
 template <typename MatrixType1, typename MatrixType2>
-inline EIGEN_DONT_INLINE bool is_approx(
-    const Eigen::MatrixBase<MatrixType1>& mat1,
-    const Eigen::MatrixBase<MatrixType2>& mat2,
-    const typename MatrixType1::Scalar& prec) {
+EIGEN_DONT_INLINE bool is_approx(const Eigen::MatrixBase<MatrixType1>& mat1,
+                                 const Eigen::MatrixBase<MatrixType2>& mat2,
+                                 const typename MatrixType1::RealScalar& prec) {
   return mat1.isApprox(mat2, prec);
 }
 
 template <typename MatrixType1, typename MatrixType2>
-inline EIGEN_DONT_INLINE bool is_approx(
-    const Eigen::MatrixBase<MatrixType1>& mat1,
-    const Eigen::MatrixBase<MatrixType2>& mat2) {
+EIGEN_DONT_INLINE bool is_approx(const Eigen::MatrixBase<MatrixType1>& mat1,
+                                 const Eigen::MatrixBase<MatrixType2>& mat2) {
   return is_approx(
       mat1, mat2,
-      Eigen::NumTraits<typename MatrixType1::Scalar>::dummy_precision());
+      Eigen::NumTraits<typename MatrixType1::RealScalar>::dummy_precision());
+}
+
+template <typename MatrixType1, typename MatrixType2>
+EIGEN_DONT_INLINE bool is_approx(
+    const Eigen::SparseMatrixBase<MatrixType1>& mat1,
+    const Eigen::SparseMatrixBase<MatrixType2>& mat2,
+    const typename MatrixType1::RealScalar& prec) {
+  return mat1.isApprox(mat2, prec);
+}
+
+template <typename MatrixType1, typename MatrixType2>
+EIGEN_DONT_INLINE bool is_approx(
+    const Eigen::SparseMatrixBase<MatrixType1>& mat1,
+    const Eigen::SparseMatrixBase<MatrixType2>& mat2) {
+  return is_approx(
+      mat1, mat2,
+      Eigen::NumTraits<typename MatrixType1::RealScalar>::dummy_precision());
 }
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_utils_scalar_is_approx_hpp__
+#endif  // ifndef __eigenpy_utils_is_approx_hpp__
diff --git a/python/main.cpp b/python/main.cpp
index 4aa48a648..6d39d1612 100644
--- a/python/main.cpp
+++ b/python/main.cpp
@@ -1,6 +1,6 @@
 /*
  * Copyright 2014-2019, CNRS
- * Copyright 2018-2023, INRIA
+ * Copyright 2018-2024, INRIA
  */
 
 #include <boost/python/scope.hpp>
@@ -17,6 +17,37 @@
 
 using namespace eigenpy;
 
+template <typename Scalar>
+void exposeIsApprox() {
+  enum { Options = 0 };
+  EIGENPY_MAKE_TYPEDEFS(Scalar, Options, s, Eigen::Dynamic, X);
+  EIGENPY_UNUSED_TYPE(VectorXs);
+  EIGENPY_UNUSED_TYPE(RowVectorXs);
+  //  typedef Eigen::SparseMatrix<Scalar, Options> SparseMatrixXs;
+  typedef typename MatrixXs::RealScalar RealScalar;
+
+  using namespace Eigen;
+  const RealScalar dummy_precision =
+      Eigen::NumTraits<RealScalar>::dummy_precision();
+
+  bp::def("is_approx",
+          (bool (*)(const Eigen::MatrixBase<MatrixXs> &,
+                    const Eigen::MatrixBase<MatrixXs> &, const RealScalar &)) &
+              is_approx,
+          (bp::arg("A"), bp::arg("B"), bp::arg("prec") = dummy_precision),
+          "Returns True if A is approximately equal to B, within the "
+          "precision determined by prec.");
+
+  //  bp::def("is_approx",
+  //          (bool (*)(const Eigen::SparseMatrixBase<SparseMatrixXs> &,
+  //                    const Eigen::SparseMatrixBase<SparseMatrixXs> &,
+  //                    const RealScalar &)) &
+  //              is_approx,
+  //          (bp::arg("A"), bp::arg("B"), bp::arg("prec") = dummy_precision),
+  //          "Returns True if A is approximately equal to B, within the "
+  //          "precision determined by prec.");
+}
+
 BOOST_PYTHON_MODULE(eigenpy_pywrap) {
   enableEigenPy();
 
@@ -46,17 +77,8 @@ BOOST_PYTHON_MODULE(eigenpy_pywrap) {
     register_symbolic_link_to_registered_type<Eigen::ComputationInfo>();
   }
 
-  {
-    using namespace Eigen;
-
-    bp::def("is_approx",
-            (bool (*)(const Eigen::MatrixBase<MatrixXd> &,
-                      const Eigen::MatrixBase<MatrixXd> &, const double &)) &
-                is_approx<MatrixXd, MatrixXd>,
-            (bp::arg("A"), bp::arg("B"), bp::arg("prec") = 1e-12),
-            "Returns True if A is approximately equal to B, within the "
-            "precision determined by prec.");
-  }
+  exposeIsApprox<double>();
+  exposeIsApprox<std::complex<double> >();
 
   exposeDecompositions();
 }
diff --git a/src/decompositions/accelerate.cpp b/src/decompositions/accelerate.cpp
new file mode 100644
index 000000000..cf0b59b2d
--- /dev/null
+++ b/src/decompositions/accelerate.cpp
@@ -0,0 +1,53 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#include "eigenpy/fwd.hpp"
+#include "eigenpy/decompositions/decompositions.hpp"
+
+#include "eigenpy/decompositions/sparse/accelerate/accelerate.hpp"
+
+namespace eigenpy {
+
+void exposeAccelerate() {
+  using namespace Eigen;
+
+  typedef Eigen::SparseMatrix<double, Eigen::ColMajor> ColMajorSparseMatrix;
+  //  typedef Eigen::SparseMatrix<double,Eigen::RowMajor> RowMajorSparseMatrix;
+
+  bp::enum_<SparseOrder_t>("SparseOrder")
+      .value("SparseOrderUser", SparseOrderUser)
+      .value("SparseOrderAMD", SparseOrderAMD)
+      .value("SparseOrderMetis", SparseOrderMetis)
+      .value("SparseOrderCOLAMD", SparseOrderCOLAMD);
+
+#define EXPOSE_ACCELERATE_DECOMPOSITION(name, doc)            \
+  AccelerateImplVisitor<name<ColMajorSparseMatrix> >::expose( \
+      EIGENPY_STRINGIZE(name), doc)
+
+  EXPOSE_ACCELERATE_DECOMPOSITION(
+      AccelerateLLT,
+      "A direct Cholesky (LLT) factorization and solver based on Accelerate.");
+  EXPOSE_ACCELERATE_DECOMPOSITION(AccelerateLDLT,
+                                  "The default Cholesky (LDLT) factorization "
+                                  "and solver based on Accelerate.");
+  EXPOSE_ACCELERATE_DECOMPOSITION(
+      AccelerateLDLTUnpivoted,
+      "A direct Cholesky-like LDL^T factorization and solver based on "
+      "Accelerate with only 1x1 pivots and no pivoting.");
+  EXPOSE_ACCELERATE_DECOMPOSITION(
+      AccelerateLDLTSBK,
+      "A direct Cholesky (LDLT) factorization and solver based on Accelerate "
+      "with Supernode Bunch-Kaufman and static pivoting.");
+  EXPOSE_ACCELERATE_DECOMPOSITION(
+      AccelerateLDLTTPP,
+      "A direct Cholesky (LDLT) factorization and solver based on Accelerate "
+      "with full threshold partial pivoting.");
+  EXPOSE_ACCELERATE_DECOMPOSITION(
+      AccelerateQR, "A QR factorization and solver based on Accelerate.");
+  EXPOSE_ACCELERATE_DECOMPOSITION(
+      AccelerateCholeskyAtA,
+      "A QR factorization and solver based on Accelerate without storing Q "
+      "(equivalent to A^TA = R^T R).");
+}
+}  // namespace eigenpy
diff --git a/src/decompositions/cholmod.cpp b/src/decompositions/cholmod.cpp
new file mode 100644
index 000000000..31237e1de
--- /dev/null
+++ b/src/decompositions/cholmod.cpp
@@ -0,0 +1,33 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#include "eigenpy/fwd.hpp"
+#include "eigenpy/decompositions/decompositions.hpp"
+
+#include "eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLLT.hpp"
+#include "eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLDLT.hpp"
+#include "eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp"
+
+namespace eigenpy {
+
+void exposeCholmod() {
+  using namespace Eigen;
+
+  typedef Eigen::SparseMatrix<double, Eigen::ColMajor> ColMajorSparseMatrix;
+  //  typedef Eigen::SparseMatrix<double,Eigen::RowMajor> RowMajorSparseMatrix;
+
+  bp::enum_<CholmodMode>("CholmodMode")
+      .value("CholmodAuto", CholmodAuto)
+      .value("CholmodSimplicialLLt", CholmodSimplicialLLt)
+      .value("CholmodSupernodalLLt", CholmodSupernodalLLt)
+      .value("CholmodLDLt", CholmodLDLt);
+
+  CholmodSimplicialLLTVisitor<ColMajorSparseMatrix>::expose(
+      "CholmodSimplicialLLT");
+  CholmodSimplicialLDLTVisitor<ColMajorSparseMatrix>::expose(
+      "CholmodSimplicialLDLT");
+  CholmodSupernodalLLTVisitor<ColMajorSparseMatrix>::expose(
+      "CholmodSupernodalLLT");
+}
+}  // namespace eigenpy
diff --git a/src/decompositions/decompositions.cpp b/src/decompositions/decompositions.cpp
index 50ebf97cc..2b58c03d0 100644
--- a/src/decompositions/decompositions.cpp
+++ b/src/decompositions/decompositions.cpp
@@ -1,5 +1,5 @@
 /*
- * Copyright 2020-2021 INRIA
+ * Copyright 2020-2024 INRIA
  */
 
 #include "eigenpy/decompositions/decompositions.hpp"
@@ -7,6 +7,9 @@
 #include "eigenpy/decompositions/EigenSolver.hpp"
 #include "eigenpy/decompositions/LDLT.hpp"
 #include "eigenpy/decompositions/LLT.hpp"
+#include "eigenpy/decompositions/PermutationMatrix.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 +18,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 +40,21 @@ void exposeDecompositions() {
         .value("ABx_lx", ABx_lx)
         .value("BAx_lx", BAx_lx);
   }
+
+  // Expose sparse decompositions
+  {
+    SimplicialLLTVisitor<ColMajorSparseMatrix>::expose("SimplicialLLT");
+    SimplicialLDLTVisitor<ColMajorSparseMatrix>::expose("SimplicialLDLT");
+  }
+
+  PermutationMatrixVisitor<Eigen::Dynamic>::expose("PermutationMatrix");
+
+#ifdef EIGENPY_WITH_CHOLMOD_SUPPORT
+  exposeCholmod();
+#endif
+
+#ifdef EIGENPY_WITH_ACCELERATE_SUPPORT
+  exposeAccelerate();
+#endif
 }
 }  // namespace eigenpy
diff --git a/src/scipy-type.cpp b/src/scipy-type.cpp
index e189e636f..ae1ae8476 100644
--- a/src/scipy-type.cpp
+++ b/src/scipy-type.cpp
@@ -28,7 +28,13 @@ const PyTypeObject* ScipyType::getScipyCSCMatrixType() {
 }
 
 ScipyType::ScipyType() {
-  sparse_module = bp::import("scipy.sparse");
+  try {
+    sparse_module = bp::import("scipy.sparse");
+  } catch (...) {
+    throw std::runtime_error(
+        "SciPy is not installed. "
+        "You can install it using the command \'pip install scipy\'.");
+  }
 
 #if PY_MAJOR_VERSION >= 3
   // TODO I don't know why this Py_INCREF is necessary.


Test patch:
diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt
index 8dc425fa9..06697b157 100644
--- a/unittest/CMakeLists.txt
+++ b/unittest/CMakeLists.txt
@@ -24,6 +24,7 @@ macro(ADD_LIB_UNIT_TEST test)
   endif(NOT BUILD_TESTING)
 endmacro(ADD_LIB_UNIT_TEST)
 
+add_dependencies(build_tests ${PYWRAP})
 add_lib_unit_test(matrix)
 if(BUILD_TESTING_SCIPY)
   find_scipy()
@@ -78,10 +79,6 @@ endif()
 add_lib_unit_test(bind_virtual_factory)
 
 add_python_unit_test("py-matrix" "unittest/python/test_matrix.py" "unittest")
-if(BUILD_TESTING_SCIPY)
-  add_python_unit_test("py-sparse-matrix"
-                       "unittest/python/test_sparse_matrix.py" "unittest")
-endif()
 
 add_python_unit_test("py-tensor" "unittest/python/test_tensor.py" "unittest")
 add_python_unit_test("py-geometry" "unittest/python/test_geometry.py"
@@ -145,3 +142,40 @@ add_python_unit_test("py-std-unique-ptr"
 
 add_python_unit_test("py-bind-virtual" "unittest/python/test_bind_virtual.py"
                      "unittest")
+
+if(BUILD_TESTING_SCIPY)
+  add_python_unit_test("py-sparse-matrix"
+                       "unittest/python/test_sparse_matrix.py" "unittest")
+
+  add_python_unit_test(
+    "py-SimplicialLLT"
+    "unittest/python/decompositions/sparse/test_SimplicialLLT.py" "python")
+  add_python_unit_test(
+    "py-SimplicialLDLT"
+    "unittest/python/decompositions/sparse/test_SimplicialLDLT.py" "python")
+
+  if(BUILD_WITH_CHOLMOD_SUPPORT)
+
+    add_python_unit_test(
+      "py-CholmodSimplicialLLT"
+      "unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLLT.py"
+      "python")
+
+    add_python_unit_test(
+      "py-CholmodSimplicialLDLT"
+      "unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLDLT.py"
+      "python")
+
+    add_python_unit_test(
+      "py-CholmodSupernodalLLT"
+      "unittest/python/decompositions/sparse/cholmod/test_CholmodSupernodalLLT.py"
+      "python")
+
+  endif(BUILD_WITH_CHOLMOD_SUPPORT)
+
+  if(BUILD_WITH_ACCELERATE_SUPPORT)
+    add_python_unit_test(
+      "py-Accelerate"
+      "unittest/python/decompositions/sparse/test_Accelerate.py" "python")
+  endif(BUILD_WITH_ACCELERATE_SUPPORT)
+endif()
diff --git a/unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLDLT.py b/unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLDLT.py
new file mode 100644
index 000000000..328518a53
--- /dev/null
+++ b/unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLDLT.py
@@ -0,0 +1,23 @@
+import eigenpy
+
+import numpy as np
+from scipy.sparse import csc_matrix
+
+dim = 100
+A = np.random.rand(dim, dim)
+A = (A + A.T) * 0.5 + np.diag(10.0 + np.random.rand(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.CholmodSimplicialLDLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+X = np.random.rand(dim, 20)
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)
diff --git a/unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLLT.py b/unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLLT.py
new file mode 100644
index 000000000..f29937705
--- /dev/null
+++ b/unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLLT.py
@@ -0,0 +1,23 @@
+import eigenpy
+
+import numpy as np
+from scipy.sparse import csc_matrix
+
+dim = 100
+A = np.random.rand(dim, dim)
+A = (A + A.T) * 0.5 + np.diag(10.0 + np.random.rand(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.CholmodSimplicialLLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+X = np.random.rand(dim, 20)
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)
diff --git a/unittest/python/decompositions/sparse/cholmod/test_CholmodSupernodalLLT.py b/unittest/python/decompositions/sparse/cholmod/test_CholmodSupernodalLLT.py
new file mode 100644
index 000000000..1fd9a7200
--- /dev/null
+++ b/unittest/python/decompositions/sparse/cholmod/test_CholmodSupernodalLLT.py
@@ -0,0 +1,23 @@
+import eigenpy
+
+import numpy as np
+from scipy.sparse import csc_matrix
+
+dim = 100
+A = np.random.rand(dim, dim)
+A = (A + A.T) * 0.5 + np.diag(10.0 + np.random.rand(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.CholmodSupernodalLLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+X = np.random.rand(dim, 20)
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)
diff --git a/unittest/python/decompositions/sparse/test_Accelerate.py b/unittest/python/decompositions/sparse/test_Accelerate.py
new file mode 100644
index 000000000..2f89a4b7d
--- /dev/null
+++ b/unittest/python/decompositions/sparse/test_Accelerate.py
@@ -0,0 +1,35 @@
+import eigenpy
+
+import numpy as np
+from scipy.sparse import csc_matrix
+
+
+def test(SolverType: type):
+    dim = 100
+    A = np.random.rand(dim, dim)
+    A = (A + A.T) * 0.5 + np.diag(10.0 + np.random.rand(dim))
+
+    A = csc_matrix(A)
+
+    llt = SolverType(A)
+
+    assert llt.info() == eigenpy.ComputationInfo.Success
+
+    X = np.random.rand(dim, 20)
+    B = A.dot(X)
+    X_est = llt.solve(B)
+    #    import pdb; pdb.set_trace()
+    assert eigenpy.is_approx(X, X_est)
+    assert eigenpy.is_approx(A.dot(X_est), B)
+
+    llt.analyzePattern(A)
+    llt.factorize(A)
+
+
+test(eigenpy.AccelerateLLT)
+test(eigenpy.AccelerateLDLT)
+test(eigenpy.AccelerateLDLTUnpivoted)
+test(eigenpy.AccelerateLDLTSBK)
+test(eigenpy.AccelerateLDLTTPP)
+test(eigenpy.AccelerateQR)
+# test(eigenpy.AccelerateCholeskyAtA) # This test is not passing. Seems there is a bug in Eigen with the support of Accelerate.
diff --git a/unittest/python/decompositions/sparse/test_SimplicialLDLT.py b/unittest/python/decompositions/sparse/test_SimplicialLDLT.py
new file mode 100644
index 000000000..b83d46a0e
--- /dev/null
+++ b/unittest/python/decompositions/sparse/test_SimplicialLDLT.py
@@ -0,0 +1,31 @@
+import eigenpy
+
+import numpy as np
+from scipy.sparse import csc_matrix
+
+dim = 100
+A = np.random.rand(dim, dim)
+A = (A + A.T) * 0.5 + np.diag(10.0 + np.random.rand(dim))
+
+A = csc_matrix(A)
+
+ldlt = eigenpy.SimplicialLDLT(A)
+
+assert ldlt.info() == eigenpy.ComputationInfo.Success
+
+L = ldlt.matrixL()
+U = ldlt.matrixU()
+D = csc_matrix(np.diag(ldlt.vectorD()))
+
+LDU = L @ D @ U
+assert eigenpy.is_approx(LDU.toarray(), A.toarray())
+
+X = np.random.rand(dim, 20)
+B = A.dot(X)
+X_est = ldlt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+ldlt.analyzePattern(A)
+ldlt.factorize(A)
+permutation = ldlt.permutationP()
diff --git a/unittest/python/decompositions/sparse/test_SimplicialLLT.py b/unittest/python/decompositions/sparse/test_SimplicialLLT.py
new file mode 100644
index 000000000..c436fea8b
--- /dev/null
+++ b/unittest/python/decompositions/sparse/test_SimplicialLLT.py
@@ -0,0 +1,30 @@
+import eigenpy
+
+import numpy as np
+from scipy.sparse import csc_matrix
+
+dim = 100
+A = np.random.rand(dim, dim)
+A = (A + A.T) * 0.5 + np.diag(10.0 + np.random.rand(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.SimplicialLLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+L = llt.matrixL()
+U = llt.matrixU()
+
+LU = L @ U
+assert eigenpy.is_approx(LU.toarray(), A.toarray())
+
+X = np.random.rand(dim, 20)
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)
+permutation = llt.permutationP()
diff --git a/unittest/user_type.cpp b/unittest/user_type.cpp
index 112dd238a..dcafea1a2 100644
--- a/unittest/user_type.cpp
+++ b/unittest/user_type.cpp
@@ -54,6 +54,9 @@ struct NumTraits<CustomType<Scalar> > {
   }
 
   static int digits10() { return std::numeric_limits<Scalar>::digits10; }
+  static int max_digits10() {
+    return std::numeric_limits<Scalar>::max_digits10;
+  }
 };
 }  // namespace Eigen