@@ -7,123 +7,144 @@ namespace nanoeigenpy {
77namespace nb = nanobind;
88
99template <typename MatrixType, typename MatrixOrVector>
10- MatrixOrVector solve (const Eigen::LDLT<MatrixType> &c, const MatrixOrVector &vec) {
11- return c.solve (vec);
10+ MatrixOrVector solve (const Eigen::LDLT<MatrixType> &c,
11+ const MatrixOrVector &vec) {
12+ return c.solve (vec);
1213}
1314
1415template <typename MatrixType>
1516void exposeLDLTSolver (nb::module_ m, const char *name) {
1617 using Solver = Eigen::LDLT<MatrixType>;
1718 using Scalar = typename MatrixType::Scalar;
1819 using VectorType = Eigen::Matrix<Scalar, -1 , 1 >;
19- auto cl = nb::class_<Solver>(m, name,
20- " Robust Cholesky decomposition of a matrix with pivoting.\n\n "
21- " Perform a robust Cholesky decomposition of a positive semidefinite or "
22- " negative semidefinite matrix $ A $ such that $ A = P^TLDL^*P $, where "
23- " P is a permutation matrix, L is lower triangular with a unit diagonal "
24- " and D is a diagonal matrix.\n\n "
25- " The decomposition uses pivoting to ensure stability, so that L will "
26- " have zeros in the bottom right rank(A) - n submatrix. Avoiding the "
27- " square root on D also stabilizes the computation." )
28-
29- .def (nb::init<>(),
30- " Default constructor." )
31- .def (nb::init<Eigen::DenseIndex>(), nb::arg (" size" ),
32- " Default constructor with memory preallocation." )
33- .def (nb::init<const MatrixType &>(), nb::arg (" matrix" ),
34- " Constructs a LLT factorization from a given matrix." )
35-
36- .def (EigenBaseVisitor ())
37-
38- .def (" isNegative" , &Solver::isNegative,
39- " Returns true if the matrix is negative (semidefinite)." )
40- .def (" isPositive" , &Solver::isPositive,
41- " Returns true if the matrix is positive (semidefinite)." )
42-
43- .def (" matrixL" ,
44- [](Solver const &c) -> MatrixType { return c.matrixL (); },
45- " Returns the lower triangular matrix L." )
46- .def (" matrixU" ,
47- [](Solver const &c) -> MatrixType { return c.matrixU (); },
48- " Returns the upper triangular matrix U." )
49- .def (" vectorD" ,
50- [](Solver const &c) -> VectorType { return c.vectorD (); },
51- " Returns the coefficients of the diagonal matrix D." )
52- .def (" matrixLDLT" , &Solver::matrixLDLT,
53- " Returns the LDLT decomposition matrix." ,
54- nb::rv_policy::reference_internal)
55-
56- .def (" transpositionsP" ,
57- [](Solver const &c) -> MatrixType { return c.transpositionsP () *
58- MatrixType::Identity (c.matrixL ().rows (), c.matrixL ().rows ()); },
59- " Returns the permutation matrix P." )
60-
61- .def (" rankUpdate" ,
62- [](Solver &c, VectorType const &w, Scalar sigma) {
63- return c.rankUpdate (w, sigma);
64- },
65- nb::arg (" w" ), nb::arg (" sigma" ))
20+ auto cl =
21+ nb::class_<Solver>(
22+ m, name,
23+ " Robust Cholesky decomposition of a matrix with pivoting.\n\n "
24+ " Perform a robust Cholesky decomposition of a positive semidefinite "
25+ " or "
26+ " negative semidefinite matrix $ A $ such that $ A = P^TLDL^*P $, "
27+ " where "
28+ " P is a permutation matrix, L is lower triangular with a unit "
29+ " diagonal "
30+ " and D is a diagonal matrix.\n\n "
31+ " The decomposition uses pivoting to ensure stability, so that L will "
32+ " have zeros in the bottom right rank(A) - n submatrix. Avoiding the "
33+ " square root on D also stabilizes the computation." )
34+
35+ .def (nb::init<>(), " Default constructor." )
36+ .def (nb::init<Eigen::DenseIndex>(), nb::arg (" size" ),
37+ " Default constructor with memory preallocation." )
38+ .def (nb::init<const MatrixType &>(), nb::arg (" matrix" ),
39+ " Constructs a LLT factorization from a given matrix." )
40+
41+ .def (EigenBaseVisitor ())
42+
43+ .def (" isNegative" , &Solver::isNegative,
44+ " Returns true if the matrix is negative (semidefinite)." )
45+ .def (" isPositive" , &Solver::isPositive,
46+ " Returns true if the matrix is positive (semidefinite)." )
47+
48+ .def (
49+ " matrixL" ,
50+ [](Solver const &c) -> MatrixType { return c.matrixL (); },
51+ " Returns the lower triangular matrix L." )
52+ .def (
53+ " matrixU" ,
54+ [](Solver const &c) -> MatrixType { return c.matrixU (); },
55+ " Returns the upper triangular matrix U." )
56+ .def (
57+ " vectorD" ,
58+ [](Solver const &c) -> VectorType { return c.vectorD (); },
59+ " Returns the coefficients of the diagonal matrix D." )
60+ .def (" matrixLDLT" , &Solver::matrixLDLT,
61+ " Returns the LDLT decomposition matrix." ,
62+ nb::rv_policy::reference_internal)
63+
64+ .def (
65+ " transpositionsP" ,
66+ [](Solver const &c) -> MatrixType {
67+ return c.transpositionsP () *
68+ MatrixType::Identity (c.matrixL ().rows (),
69+ c.matrixL ().rows ());
70+ },
71+ " Returns the permutation matrix P." )
72+
73+ .def (
74+ " rankUpdate" ,
75+ [](Solver &c, VectorType const &w, Scalar sigma) {
76+ return c.rankUpdate (w, sigma);
77+ },
78+ nb::arg (" w" ), nb::arg (" sigma" ))
6679
6780#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
68- .def (" adjoint" , &Solver::adjoint,
69- " Returns the adjoint, that is, a reference to the decomposition "
70- " itself as if the underlying matrix is self-adjoint." ,
71- nb::rv_policy::reference)
81+ .def (" adjoint" , &Solver::adjoint,
82+ " Returns the adjoint, that is, a reference to the decomposition "
83+ " itself as if the underlying matrix is self-adjoint." ,
84+ nb::rv_policy::reference)
7285#endif
7386
74- .def (" compute " ,
75- [](Solver &c, VectorType const &matrix) {
76- return c. compute ( matrix);
77- },
78- nb::arg ( " matrix " ) ,
79- " Computes the LDLT of given matrix." ,
80- nb::rv_policy::reference)
81- .def (" info" , &Solver::info,
82- " NumericalIssue if the input contains INF or NaN values or "
83- " overflow occured. Returns Success otherwise." )
87+ .def (
88+ " compute " ,
89+ [](Solver &c, VectorType const & matrix) {
90+ return c. compute (matrix);
91+ } ,
92+ nb::arg ( " matrix " ), " Computes the LDLT of given matrix." ,
93+ nb::rv_policy::reference)
94+ .def (" info" , &Solver::info,
95+ " NumericalIssue if the input contains INF or NaN values or "
96+ " overflow occured. Returns Success otherwise." )
8497
8598#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
86- .def (" rcond" , &Solver::rcond,
87- " Returns an estimate of the reciprocal condition number of the "
88- " matrix." )
99+ .def (" rcond" , &Solver::rcond,
100+ " Returns an estimate of the reciprocal condition number of the "
101+ " matrix." )
89102#endif
90103
91- .def (" reconstructedMatrix" , &Solver::reconstructedMatrix,
92- " Returns the matrix represented by the decomposition, i.e., it "
93- " returns the product: L L^*. This function is provided for debug "
94- " purpose." )
95-
96- .def (" solve" ,
97- [](Solver const &c, VectorType const &b) -> VectorType { return solve (c, b); },
98- nb::arg (" b" ),
99- " Returns the solution x of A x = b using the current "
100- " decomposition of A." )
101- .def (" solve" ,
102- [](Solver const &c, MatrixType const &B) -> MatrixType { return solve (c, B); },
103- nb::arg (" B" ),
104- " Returns the solution X of A X = B using the current "
105- " decomposition of A where B is a right hand side matrix." )
106-
107- .def (" setZero" , &Solver::setZero,
108- " Clear any existing decomposition." )
109-
110- .def (" id" ,
111- [](Solver const &c) -> int64_t { return reinterpret_cast <int64_t >(&c); },
112- " Returns the unique identity of an object.\n "
113- " For objects held in C++, it corresponds to its memory address." );
114-
104+ .def (
105+ " reconstructedMatrix" , &Solver::reconstructedMatrix,
106+ " Returns the matrix represented by the decomposition, i.e., it "
107+ " returns the product: L L^*. This function is provided for debug "
108+ " purpose." )
109+
110+ .def (
111+ " solve" ,
112+ [](Solver const &c, VectorType const &b) -> VectorType {
113+ return solve (c, b);
114+ },
115+ nb::arg (" b" ),
116+ " Returns the solution x of A x = b using the current "
117+ " decomposition of A." )
118+ .def (
119+ " solve" ,
120+ [](Solver const &c, MatrixType const &B) -> MatrixType {
121+ return solve (c, B);
122+ },
123+ nb::arg (" B" ),
124+ " Returns the solution X of A X = B using the current "
125+ " decomposition of A where B is a right hand side matrix." )
126+
127+ .def (" setZero" , &Solver::setZero, " Clear any existing decomposition." )
128+
129+ .def (
130+ " id" ,
131+ [](Solver const &c) -> int64_t {
132+ return reinterpret_cast <int64_t >(&c);
133+ },
134+ " Returns the unique identity of an object.\n "
135+ " For objects held in C++, it corresponds to its memory address." );
115136}
116137
117138} // namespace nanoeigenpy
118139
119-
120140// TODO
121141
122- // Tests that were not done in eigenpy that we could add in nanoeigenpy: (+ those cited in llt.hpp)
123- // setZero
142+ // Tests that were not done in eigenpy that we could add in nanoeigenpy: (+
143+ // those cited in llt.hpp) setZero
124144
125145// Expose supplementary content:
126146// setZero in LLT decomp too ? (not done in eigenpy)
127147
128148// Questions about eigenpy itself:
129- // Relevant to have the reconstructedMatrix method ? (knowing that we are on LDLT, not LLT)
149+ // Relevant to have the reconstructedMatrix method ? (knowing that we are on
150+ // LDLT, not LLT)
0 commit comments