Skip to content

Commit a4bfa94

Browse files
jcarpentjcarpent
authored
Merge pull request #438 from jcarpent/topic/sparse
Expose Simplicial{LLT,LDLT}
2 parents 6a0ebc5 + 37721b5 commit a4bfa94

34 files changed

+1283
-33
lines changed

‎CHANGELOG.md

+1
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
1212
- Support for `std::variant` types with `VariantConverter` ([#431](https://github.com/stack-of-tasks/eigenpy/pull/431))
1313
- 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))
1414
- Support for `std::unique_ptr` as an internal reference with `ReturnInternalStdUniquePtr` ([#433](https://github.com/stack-of-tasks/eigenpy/pull/433))
15+
- 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))
1516

1617
### Fixed
1718
- Fix the issue of missing exposition of Eigen types with __int64 scalar type ([#426](https://github.com/stack-of-tasks/eigenpy/pull/426))

‎CMakeLists.txt

+92
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,15 @@ include("${JRL_CMAKE_MODULES}/apple.cmake")
5656
option(GENERATE_PYTHON_STUBS
5757
"Generate the Python stubs associated to the Python library" OFF)
5858

59+
option(BUILD_WITH_CHOLMOD_SUPPORT "Build EigenPy with the Cholmod support" OFF)
60+
61+
if(APPLE)
62+
option(BUILD_WITH_ACCELERATE_SUPPORT
63+
"Build EigenPy with the Accelerate support" OFF)
64+
else(APPLE)
65+
set(BUILD_WITH_ACCELERATE_SUPPORT FALSE)
66+
endif(APPLE)
67+
5968
string(REPLACE "-pedantic" "" CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS})
6069

6170
# If needed, fix CMake policy for APPLE systems
@@ -99,6 +108,35 @@ export_boost_default_options()
99108
find_package(Boost REQUIRED)
100109
search_for_boost_python(REQUIRED)
101110

111+
if(BUILD_WITH_CHOLMOD_SUPPORT)
112+
set(CMAKE_MODULE_PATH ${JRL_CMAKE_MODULES}/find-external/CHOLMOD
113+
${CMAKE_MODULE_PATH})
114+
add_project_dependency(CHOLMOD REQUIRED FIND_EXTERNAL "CHOLMOD")
115+
message(
116+
STATUS
117+
"Build with CHOLDOD support (LGPL). See CHOLMOD/Doc/License.txt for further details."
118+
)
119+
add_definitions(-DEIGENPY_WITH_CHOLMOD_SUPPORT)
120+
endif(BUILD_WITH_CHOLMOD_SUPPORT)
121+
122+
if(BUILD_WITH_ACCELERATE_SUPPORT)
123+
if(NOT ${Eigen3_VERSION} VERSION_GREATER_EQUAL "3.4.90")
124+
message(
125+
FATAL_ERROR
126+
"Your version of Eigen is too low. Should be at least 3.4.90. Current version is ${Eigen3_VERSION}."
127+
)
128+
endif()
129+
130+
set(CMAKE_MODULE_PATH ${JRL_CMAKE_MODULES}/find-external/Accelerate
131+
${CMAKE_MODULE_PATH})
132+
find_package(
133+
Accelerate REQUIRED # FIND_EXTERNAL "Accelerate" # We don't export yet as
134+
# there might be an issue on AMR64 platforms
135+
)
136+
message(STATUS "Build with Accelerate support framework.")
137+
add_definitions(-DEIGENPY_WITH_ACCELERATE_SUPPORT)
138+
endif(BUILD_WITH_ACCELERATE_SUPPORT)
139+
102140
# ----------------------------------------------------
103141
# --- INCLUDE ----------------------------------------
104142
# ----------------------------------------------------
@@ -117,9 +155,39 @@ set(${PROJECT_NAME}_SOLVERS_HEADERS
117155
include/eigenpy/solvers/BasicPreconditioners.hpp
118156
include/eigenpy/solvers/BFGSPreconditioners.hpp)
119157

158+
set(${PROJECT_NAME}_EIGEN_HEADERS include/eigenpy/eigen/EigenBase.hpp)
159+
160+
set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_CHOLMOD_HEADERS
161+
include/eigenpy/decompositions/sparse/cholmod/CholmodBase.hpp
162+
include/eigenpy/decompositions/sparse/cholmod/CholmodDecomposition.hpp
163+
include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLDLT.hpp
164+
include/eigenpy/decompositions/sparse/cholmod/CholmodSimplicialLLT.hpp
165+
include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp)
166+
167+
set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_ACCELERATE_HEADERS
168+
include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp)
169+
170+
set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
171+
include/eigenpy/decompositions/sparse/LLT.hpp
172+
include/eigenpy/decompositions/sparse/LDLT.hpp
173+
include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
174+
include/eigenpy/decompositions/sparse/SparseSolverBase.hpp)
175+
176+
if(BUILD_WITH_CHOLMOD_SUPPORT)
177+
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
178+
${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_CHOLMOD_HEADERS})
179+
endif(BUILD_WITH_CHOLMOD_SUPPORT)
180+
181+
if(BUILD_WITH_ACCELERATE_SUPPORT)
182+
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
183+
${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_ACCELERATE_HEADERS})
184+
endif(BUILD_WITH_ACCELERATE_SUPPORT)
185+
120186
set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
187+
${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS}
121188
include/eigenpy/decompositions/decompositions.hpp
122189
include/eigenpy/decompositions/EigenSolver.hpp
190+
include/eigenpy/decompositions/PermutationMatrix.hpp
123191
include/eigenpy/decompositions/LDLT.hpp
124192
include/eigenpy/decompositions/LLT.hpp
125193
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
@@ -128,6 +196,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
128196
set(${PROJECT_NAME}_HEADERS
129197
${${PROJECT_NAME}_UTILS_HEADERS}
130198
${${PROJECT_NAME}_SOLVERS_HEADERS}
199+
${${PROJECT_NAME}_EIGEN_HEADERS}
131200
${${PROJECT_NAME}_DECOMPOSITIONS_HEADERS}
132201
include/eigenpy/alignment.hpp
133202
include/eigenpy/computation-info.hpp
@@ -188,6 +257,16 @@ set(${PROJECT_NAME}_SOLVERS_SOURCES src/solvers/preconditioners.cpp
188257
set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
189258
src/decompositions/decompositions.cpp)
190259

260+
if(BUILD_WITH_CHOLMOD_SUPPORT)
261+
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
262+
src/decompositions/cholmod.cpp)
263+
endif(BUILD_WITH_CHOLMOD_SUPPORT)
264+
265+
if(BUILD_WITH_ACCELERATE_SUPPORT)
266+
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
267+
src/decompositions/accelerate.cpp)
268+
endif(BUILD_WITH_ACCELERATE_SUPPORT)
269+
191270
set(${PROJECT_NAME}_SOURCES
192271
${${PROJECT_NAME}_SOLVERS_SOURCES}
193272
${${PROJECT_NAME}_DECOMPOSITIONS_SOURCES}
@@ -239,6 +318,19 @@ modernize_target_link_libraries(
239318
${NUMPY_INCLUDE_DIRS}
240319
${PYTHON_INCLUDE_DIR})
241320

321+
# Links against CholMod
322+
if(BUILD_WITH_CHOLMOD_SUPPORT)
323+
modernize_target_link_libraries(${PROJECT_NAME} SCOPE PUBLIC TARGETS
324+
CHOLMOD::CHOLMOD)
325+
endif(BUILD_WITH_CHOLMOD_SUPPORT)
326+
327+
# Links against accelerate
328+
if(BUILD_WITH_ACCELERATE_SUPPORT)
329+
# modernize_target_link_libraries(${PROJECT_NAME} SCOPE PUBLIC TARGETS
330+
# Accelerate)
331+
target_link_libraries(${PROJECT_NAME} PRIVATE "-framework accelerate")
332+
endif(BUILD_WITH_ACCELERATE_SUPPORT)
333+
242334
if(SUFFIX_SO_VERSION)
243335
set_target_properties(${PROJECT_NAME} PROPERTIES SOVERSION ${PROJECT_VERSION})
244336
endif(SUFFIX_SO_VERSION)

‎cmake

‎include/eigenpy/decompositions/LDLT.hpp

+3
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010

1111
#include "eigenpy/eigenpy.hpp"
1212
#include "eigenpy/utils/scalar-name.hpp"
13+
#include "eigenpy/eigen/EigenBase.hpp"
1314

1415
namespace eigenpy {
1516

@@ -36,6 +37,8 @@ struct LDLTSolverVisitor
3637
bp::args("self", "matrix"),
3738
"Constructs a LDLT factorization from a given matrix."))
3839

40+
.def(EigenBaseVisitor<Solver>())
41+
3942
.def("isNegative", &Solver::isNegative, bp::arg("self"),
4043
"Returns true if the matrix is negative (semidefinite).")
4144
.def("isPositive", &Solver::isPositive, bp::arg("self"),

‎include/eigenpy/decompositions/LLT.hpp

+3
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010

1111
#include "eigenpy/eigenpy.hpp"
1212
#include "eigenpy/utils/scalar-name.hpp"
13+
#include "eigenpy/eigen/EigenBase.hpp"
1314

1415
namespace eigenpy {
1516

@@ -36,6 +37,8 @@ struct LLTSolverVisitor
3637
bp::args("self", "matrix"),
3738
"Constructs a LLT factorization from a given matrix."))
3839

40+
.def(EigenBaseVisitor<Solver>())
41+
3942
.def("matrixL", &matrixL, bp::arg("self"),
4043
"Returns the lower triangular matrix L.")
4144
.def("matrixU", &matrixU, bp::arg("self"),

‎include/eigenpy/decompositions/PermutationMatrix.hpp

+106
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decomposition_permutation_matrix_hpp__
6+
#define __eigenpy_decomposition_permutation_matrix_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
#include "eigenpy/eigen/EigenBase.hpp"
10+
11+
namespace eigenpy {
12+
13+
template <int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime,
14+
typename StorageIndex_ = int>
15+
struct PermutationMatrixVisitor
16+
: public boost::python::def_visitor<PermutationMatrixVisitor<
17+
SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> > {
18+
typedef StorageIndex_ StorageIndex;
19+
typedef Eigen::PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime,
20+
StorageIndex>
21+
PermutationMatrix;
22+
typedef typename PermutationMatrix::DenseMatrixType DenseMatrixType;
23+
typedef PermutationMatrix Self;
24+
typedef Eigen::Matrix<StorageIndex, SizeAtCompileTime, 1, 0,
25+
MaxSizeAtCompileTime, 1>
26+
VectorIndex;
27+
28+
template <class PyClass>
29+
void visit(PyClass &cl) const {
30+
cl.def(bp::init<const Eigen::DenseIndex>(bp::args("self", "size"),
31+
"Default constructor"))
32+
.def(bp::init<VectorIndex>(
33+
bp::args("self", "indices"),
34+
"The indices array has the meaning that the permutations sends "
35+
"each integer i to indices[i].\n"
36+
"It is your responsibility to check that the indices array that "
37+
"you passes actually describes a permutation, i.e., each value "
38+
"between 0 and n-1 occurs exactly once, where n is the array's "
39+
"size."))
40+
41+
.def(
42+
"indices",
43+
+[](const PermutationMatrix &self) {
44+
return VectorIndex(self.indices());
45+
},
46+
bp::arg("self"), "The stored array representing the permutation.")
47+
48+
.def("applyTranspositionOnTheLeft",
49+
&PermutationMatrix::applyTranspositionOnTheLeft,
50+
bp::args("self", "i", "j"),
51+
"Multiplies self by the transposition (ij) on the left.",
52+
bp::return_self<>())
53+
.def("applyTranspositionOnTheRight",
54+
&PermutationMatrix::applyTranspositionOnTheRight,
55+
bp::args("self", "i", "j"),
56+
"Multiplies self by the transposition (ij) on the right.",
57+
bp::return_self<>())
58+
59+
.def("setIdentity",
60+
(void(PermutationMatrix::*)()) & PermutationMatrix::setIdentity,
61+
bp::arg("self"),
62+
"Sets self to be the identity permutation matrix.")
63+
.def("setIdentity",
64+
(void(PermutationMatrix::*)(Eigen::DenseIndex)) &
65+
PermutationMatrix::setIdentity,
66+
bp::args("self", "size"),
67+
"Sets self to be the identity permutation matrix of given size.")
68+
69+
.def("toDenseMatrix", &PermutationMatrix::toDenseMatrix,
70+
bp::arg("self"),
71+
"Returns a numpy array object initialized from this permutation "
72+
"matrix.")
73+
74+
.def(
75+
"transpose",
76+
+[](const PermutationMatrix &self) -> PermutationMatrix {
77+
return self.transpose();
78+
},
79+
bp::arg("self"), "Returns the tranpose permutation matrix.")
80+
.def(
81+
"inverse",
82+
+[](const PermutationMatrix &self) -> PermutationMatrix {
83+
return self.inverse();
84+
},
85+
bp::arg("self"), "Returns the inverse permutation matrix.")
86+
87+
.def("resize", &PermutationMatrix::resize, bp::args("self", "size"),
88+
"Resizes to given size.")
89+
90+
.def(bp::self * bp::self)
91+
.def(EigenBaseVisitor<Self>());
92+
}
93+
94+
static void expose(const std::string &name) {
95+
bp::class_<PermutationMatrix>(name.c_str(),
96+
"Permutation matrix.\n"
97+
"This class represents a permutation matrix, "
98+
"internally stored as a vector of integers.",
99+
bp::no_init)
100+
.def(PermutationMatrixVisitor());
101+
}
102+
};
103+
104+
} // namespace eigenpy
105+
106+
#endif // ifndef __eigenpy_decomposition_permutation_matrix_hpp__

‎include/eigenpy/decompositions/decompositions.hpp

+9
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,15 @@
99

1010
namespace eigenpy {
1111
void EIGENPY_DLLAPI exposeDecompositions();
12+
13+
#ifdef EIGENPY_WITH_CHOLMOD_SUPPORT
14+
void EIGENPY_DLLAPI exposeCholmod();
15+
#endif
16+
17+
#ifdef EIGENPY_WITH_ACCELERATE_SUPPORT
18+
void EIGENPY_DLLAPI exposeAccelerate();
19+
#endif
20+
1221
} // namespace eigenpy
1322

1423
#endif // define __eigenpy_decompositions_decompositions_hpp__

‎include/eigenpy/decompositions/sparse/LDLT.hpp

+71
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decomposition_sparse_ldlt_hpp__
6+
#define __eigenpy_decomposition_sparse_ldlt_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
10+
#include "eigenpy/utils/scalar-name.hpp"
11+
12+
namespace eigenpy {
13+
14+
template <typename _MatrixType, int _UpLo = Eigen::Lower,
15+
typename _Ordering =
16+
Eigen::AMDOrdering<typename _MatrixType::StorageIndex> >
17+
struct SimplicialLDLTVisitor
18+
: public boost::python::def_visitor<
19+
SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering> > {
20+
typedef SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor;
21+
typedef _MatrixType MatrixType;
22+
23+
typedef Eigen::SimplicialLDLT<MatrixType> Solver;
24+
typedef typename MatrixType::Scalar Scalar;
25+
typedef typename MatrixType::RealScalar RealScalar;
26+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
27+
DenseVectorXs;
28+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
29+
MatrixType::Options>
30+
DenseMatrixXs;
31+
32+
template <class PyClass>
33+
void visit(PyClass &cl) const {
34+
cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
35+
.def(bp::init<MatrixType>(bp::args("self", "matrix"),
36+
"Constructs and performs the LDLT "
37+
"factorization from a given matrix."))
38+
39+
.def("vectorD", &vectorD, bp::arg("self"),
40+
"Returns the diagonal vector D.")
41+
.def(SimplicialCholeskyVisitor<Solver>());
42+
}
43+
44+
static void expose() {
45+
static const std::string classname =
46+
"SimplicialLDLT_" + scalar_name<Scalar>::shortname();
47+
expose(classname);
48+
}
49+
50+
static void expose(const std::string &name) {
51+
bp::class_<Solver, boost::noncopyable>(
52+
name.c_str(),
53+
"A direct sparse LDLT Cholesky factorizations.\n\n"
54+
"This class provides a LDL^T Cholesky factorizations of sparse "
55+
"matrices that are selfadjoint and positive definite."
56+
"The factorization allows for solving A.X = B where X and B can be "
57+
"either dense or sparse.\n\n"
58+
"In order to reduce the fill-in, a symmetric permutation P is applied "
59+
"prior to the factorization such that the factorized matrix is P A "
60+
"P^-1.",
61+
bp::no_init)
62+
.def(SimplicialLDLTVisitor());
63+
}
64+
65+
private:
66+
static DenseVectorXs vectorD(const Solver &self) { return self.vectorD(); }
67+
};
68+
69+
} // namespace eigenpy
70+
71+
#endif // ifndef __eigenpy_decomposition_sparse_ldlt_hpp__

0 commit comments

Comments
 (0)