Skip to content

Commit f5da84e

Browse files
Lucas-HaubertLucas-Haubert
committed
Add IncompleteLUT
1 parent 9b97fd9 commit f5da84e

File tree

7 files changed

+239
-67
lines changed

7 files changed

+239
-67
lines changed

‎CMakeLists.txt‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -214,7 +214,8 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
214214
include/eigenpy/decompositions/sparse/LU.hpp
215215
include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
216216
include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
217-
include/eigenpy/decompositions/sparse/IncompleteCholesky.hpp)
217+
include/eigenpy/decompositions/sparse/IncompleteCholesky.hpp
218+
include/eigenpy/decompositions/sparse/IncompleteLUT.hpp)
218219

219220
if(BUILD_WITH_CHOLMOD_SUPPORT)
220221
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
@@ -342,7 +343,8 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
342343
src/decompositions/permutation-matrix.cpp
343344
src/decompositions/simplicial-llt-solver.cpp
344345
src/decompositions/simplicial-ldlt-solver.cpp
345-
src/decompositions/incomplete-cholesky.cpp)
346+
src/decompositions/incomplete-cholesky.cpp
347+
src/decompositions/incomplete-lut.cpp)
346348

347349
if(BUILD_WITH_CHOLMOD_SUPPORT)
348350
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SOURCES

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

Lines changed: 59 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -2,19 +2,17 @@
22
* Copyright 2024 INRIA
33
*/
44

5-
#ifndef __eigenpy_decompositions_sparse_ldlt_hpp__
6-
#define __eigenpy_decompositions_sparse_ldlt_hpp__
5+
#ifndef __eigenpy_decompositions_sparse_incomplete_cholesky_hpp__
6+
#define __eigenpy_decompositions_sparse_incomplete_cholesky_hpp__
77

88
#include "eigenpy/eigenpy.hpp"
99
#include "eigenpy/utils/scalar-name.hpp"
1010

1111
namespace eigenpy {
1212

1313
template <typename _MatrixType>
14-
struct IncompleteCholeskyVisitor
15-
: public boost::python::def_visitor<
16-
IncompleteCholeskyVisitor<_MatrixType>> {
17-
14+
struct IncompleteCholeskyVisitor : public boost::python::def_visitor<
15+
IncompleteCholeskyVisitor<_MatrixType>> {
1816
typedef _MatrixType MatrixType;
1917
typedef typename MatrixType::Scalar Scalar;
2018
typedef typename MatrixType::RealScalar RealScalar;
@@ -27,82 +25,81 @@ struct IncompleteCholeskyVisitor
2725

2826
typedef Eigen::SparseMatrix<Scalar, Eigen::ColMajor> FactorType;
2927
typedef Eigen::Matrix<RealScalar, Eigen::Dynamic, 1> VectorRx;
30-
typedef Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> PermutationType;
28+
typedef Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>
29+
PermutationType;
3130

3231
template <class PyClass>
33-
void visit(PyClass &cl) const {
32+
void visit(PyClass& cl) const {
3433
cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
3534
.def(bp::init<MatrixType>(bp::args("self", "matrix"),
3635
"Constructs and performs the LDLT "
3736
"factorization from a given matrix."))
3837

39-
.def("rows", &Solver::rows, bp::arg("self"),
38+
.def("rows", &Solver::rows, bp::arg("self"),
4039
"Returns the number of rows of the matrix.")
41-
.def("cols", &Solver::cols, bp::arg("self"),
40+
.def("cols", &Solver::cols, bp::arg("self"),
4241
"Returns the number of cols of the matrix.")
4342

4443
.def("info", &Solver::info, bp::arg("self"),
45-
"Reports whether previous computation was successful.")
44+
"Reports whether previous computation was successful.")
4645

47-
.def("setInitialShift", &Solver::setInitialShift,
48-
bp::args("self", "shift"),
49-
"Set the initial shift parameter.")
46+
.def("setInitialShift", &Solver::setInitialShift,
47+
bp::args("self", "shift"), "Set the initial shift parameter.")
5048

5149
.def(
52-
"analyzePattern",
53-
+[](Solver& self, const MatrixType& amat) {
54-
self.analyzePattern(amat);
55-
},
56-
bp::arg("matrix"))
50+
"analyzePattern",
51+
+[](Solver& self, const MatrixType& amat) {
52+
self.analyzePattern(amat);
53+
},
54+
bp::arg("matrix"))
5755
.def(
58-
"factorize",
59-
+[](Solver& self, const MatrixType& amat) { self.factorize(amat); },
60-
bp::arg("matrix"))
56+
"factorize",
57+
+[](Solver& self, const MatrixType& amat) { self.factorize(amat); },
58+
bp::arg("matrix"))
6159
.def(
62-
"compute",
63-
+[](Solver& self, const MatrixType& amat) { self.compute(amat); },
64-
bp::arg("matrix"))
60+
"compute",
61+
+[](Solver& self, const MatrixType& amat) { self.compute(amat); },
62+
bp::arg("matrix"))
6563

6664
.def(
67-
"matrixL",
68-
+[](Solver& self) -> const FactorType& { return self.matrixL(); },
69-
bp::return_value_policy<bp::copy_const_reference>())
65+
"matrixL",
66+
+[](Solver& self) -> const FactorType& { return self.matrixL(); },
67+
bp::return_value_policy<bp::copy_const_reference>())
7068
.def(
71-
"scalingS",
72-
+[](Solver& self) -> const VectorRx& { return self.scalingS(); },
73-
bp::return_value_policy<bp::copy_const_reference>())
69+
"scalingS",
70+
+[](Solver& self) -> const VectorRx& { return self.scalingS(); },
71+
bp::return_value_policy<bp::copy_const_reference>())
7472
.def(
75-
"permutationP",
76-
+[](Solver& self) -> const PermutationType& {
77-
return self.permutationP();
78-
},
79-
bp::return_value_policy<bp::copy_const_reference>())
73+
"permutationP",
74+
+[](Solver& self) -> const PermutationType& {
75+
return self.permutationP();
76+
},
77+
bp::return_value_policy<bp::copy_const_reference>())
8078

8179
.def(
82-
"solve",
83-
+[](Solver const& self, const Eigen::Ref<DenseVectorXs const>& b)
84-
-> DenseVectorXs { return self.solve(b); },
85-
bp::arg("b"),
86-
"Returns the solution x of A x = b using the current decomposition "
87-
"of A, where b is a right hand side vector.")
80+
"solve",
81+
+[](Solver const& self, const Eigen::Ref<DenseVectorXs const>& b)
82+
-> DenseVectorXs { return self.solve(b); },
83+
bp::arg("b"),
84+
"Returns the solution x of A x = b using the current decomposition "
85+
"of A, where b is a right hand side vector.")
8886
.def(
89-
"solve",
90-
+[](Solver const& self, const Eigen::Ref<DenseMatrixXs const>& B)
91-
-> DenseMatrixXs { return self.solve(B); },
92-
bp::arg("b"),
93-
"Returns the solution X of A X = B using the current decomposition "
94-
"of A where B is a right hand side matrix.")
87+
"solve",
88+
+[](Solver const& self, const Eigen::Ref<DenseMatrixXs const>& B)
89+
-> DenseMatrixXs { return self.solve(B); },
90+
bp::arg("b"),
91+
"Returns the solution X of A X = B using the current decomposition "
92+
"of A where B is a right hand side matrix.")
9593
.def(
96-
"solve",
97-
+[](Solver const& self, const MatrixType& B) -> MatrixType {
98-
DenseMatrixXs B_dense = DenseMatrixXs(B);
99-
DenseMatrixXs X_dense = self.solve(B_dense);
100-
return MatrixType(X_dense.sparseView());
101-
},
102-
bp::arg("b"),
103-
"Returns the solution X of A X = B using the current decomposition "
104-
"of A where B is a right hand side matrix.")
105-
;
94+
"solve",
95+
+[](Solver const& self, const MatrixType& B) -> MatrixType {
96+
DenseMatrixXs B_dense = DenseMatrixXs(B);
97+
DenseMatrixXs X_dense = self.solve(B_dense);
98+
return MatrixType(X_dense.sparseView());
99+
},
100+
bp::arg("b"),
101+
"Returns the solution X of A X = B using the current decomposition "
102+
"of A where B is a right hand side matrix.");
106103
}
107104

108105
static void expose() {
@@ -111,16 +108,14 @@ struct IncompleteCholeskyVisitor
111108
expose(classname);
112109
}
113110

114-
static void expose(const std::string &name) {
115-
bp::class_<Solver, boost::noncopyable>(
116-
name.c_str(),
117-
"Incomplete Cholesky.",
118-
bp::no_init)
111+
static void expose(const std::string& name) {
112+
bp::class_<Solver, boost::noncopyable>(name.c_str(), "Incomplete Cholesky.",
113+
bp::no_init)
119114
.def(IncompleteCholeskyVisitor())
120115
.def(IdVisitor<Solver>());
121116
}
122117
};
123118

124119
} // namespace eigenpy
125120

126-
#endif // ifndef __eigenpy_decompositions_sparse_ldlt_hpp__
121+
#endif // ifndef __eigenpy_decompositions_sparse_incomplete_cholesky_hpp__

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

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_sparse_incomplete_lut_hpp__
6+
#define __eigenpy_decompositions_sparse_incomplete_lut_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
#include "eigenpy/utils/scalar-name.hpp"
10+
11+
namespace eigenpy {
12+
13+
template <typename _MatrixType>
14+
struct IncompleteLUTVisitor
15+
: public boost::python::def_visitor<
16+
IncompleteLUTVisitor<_MatrixType>> {
17+
18+
typedef _MatrixType MatrixType;
19+
typedef typename MatrixType::Scalar Scalar;
20+
typedef typename MatrixType::RealScalar RealScalar;
21+
typedef Eigen::IncompleteLUT<Scalar> Solver;
22+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
23+
DenseVectorXs;
24+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
25+
MatrixType::Options>
26+
DenseMatrixXs;
27+
28+
template <class PyClass>
29+
void visit(PyClass &cl) const {
30+
cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
31+
.def(bp::init<MatrixType>(bp::args("self", "matrix"),
32+
"Constructs and performs the LDLT "
33+
"factorization from a given matrix."))
34+
.def(bp::init<const MatrixType&, RealScalar, int>(
35+
(bp::arg("matrix"),
36+
bp::arg("droptol") = Eigen::NumTraits<Scalar>::dummy_precision(),
37+
bp::arg("fillfactor") = 10),
38+
"Constructs an incomplete LU factorization from a given matrix."))
39+
40+
.def("rows", &Solver::rows, bp::arg("self"),
41+
"Returns the number of rows of the matrix.")
42+
.def("cols", &Solver::cols, bp::arg("self"),
43+
"Returns the number of cols of the matrix.")
44+
45+
.def("info", &Solver::info, bp::arg("self"),
46+
"Reports whether previous computation was successful.")
47+
48+
.def(
49+
"analyzePattern",
50+
+[](Solver& self, const MatrixType& amat) {
51+
self.analyzePattern(amat);
52+
},
53+
bp::arg("matrix"))
54+
.def(
55+
"factorize",
56+
+[](Solver& self, const MatrixType& amat) { self.factorize(amat); },
57+
bp::arg("matrix"))
58+
.def(
59+
"compute",
60+
+[](Solver& self, const MatrixType& amat) { self.compute(amat); },
61+
bp::arg("matrix"))
62+
63+
.def("setDroptol", &Solver::setDroptol, bp::arg("self"))
64+
.def("setFillfactor", &Solver::setFillfactor, bp::arg("self"))
65+
66+
.def(
67+
"solve",
68+
+[](Solver const& self, const Eigen::Ref<DenseVectorXs const>& b)
69+
-> DenseVectorXs { return self.solve(b); },
70+
bp::arg("b"),
71+
"Returns the solution x of A x = b using the current decomposition "
72+
"of A, where b is a right hand side vector.")
73+
.def(
74+
"solve",
75+
+[](Solver const& self, const Eigen::Ref<DenseMatrixXs const>& B)
76+
-> DenseMatrixXs { return self.solve(B); },
77+
bp::arg("b"),
78+
"Returns the solution X of A X = B using the current decomposition "
79+
"of A where B is a right hand side matrix.")
80+
.def(
81+
"solve",
82+
+[](Solver const& self, const MatrixType& B) -> MatrixType {
83+
DenseMatrixXs B_dense = DenseMatrixXs(B);
84+
DenseMatrixXs X_dense = self.solve(B_dense);
85+
return MatrixType(X_dense.sparseView());
86+
},
87+
bp::arg("b"),
88+
"Returns the solution X of A X = B using the current decomposition "
89+
"of A where B is a right hand side matrix.")
90+
;
91+
}
92+
93+
static void expose() {
94+
static const std::string classname =
95+
"IncompleteLUT_" + scalar_name<Scalar>::shortname();
96+
expose(classname);
97+
}
98+
99+
static void expose(const std::string &name) {
100+
bp::class_<Solver, boost::noncopyable>(
101+
name.c_str(),
102+
"Incomplete LU factorization with dual-threshold strategy.",
103+
bp::no_init)
104+
.def(IncompleteLUTVisitor())
105+
.def(IdVisitor<Solver>());
106+
}
107+
};
108+
109+
} // namespace eigenpy
110+
111+
#endif // ifndef __eigenpy_decompositions_sparse_incomplete_lut_hpp__

‎src/decompositions/decompositions.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ void exposeMINRESSolver();
2727
void exposeSimplicialLLTSolver();
2828
void exposeSimplicialLDLTSolver();
2929
void exposeIncompleteCholesky();
30+
void exposeIncompleteLUT();
3031
void exposeSparseLUSolver();
3132
void exposeSparseQRSolver();
3233
void exposePermutationMatrix();
@@ -74,6 +75,7 @@ void exposeDecompositions() {
7475
exposeSparseLUSolver();
7576
exposeSparseQRSolver();
7677
exposeIncompleteCholesky();
78+
exposeIncompleteLUT();
7779

7880
exposePermutationMatrix();
7981

‎src/decompositions/incomplete-lut.cpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#include "eigenpy/decompositions/sparse/IncompleteLUT.hpp"
6+
7+
namespace eigenpy {
8+
void exposeIncompleteLUT() {
9+
using namespace Eigen;
10+
typedef SparseMatrix<double, ColMajor> ColMajorSparseMatrix;
11+
IncompleteLUTVisitor<ColMajorSparseMatrix>::expose("IncompleteLUT");
12+
}
13+
} // namespace eigenpy

‎unittest/python/decompositions/sparse/test_IncompleteCholesky.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
import numpy as np
22
from scipy.sparse import csc_matrix
3-
import eigenpy
43

4+
import eigenpy
55

66
dim = 100
77
rng = np.random.default_rng()

‎unittest/python/decompositions/sparse/test_IncompleteLUT.py‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
import numpy as np
2+
from scipy.sparse import csc_matrix
3+
import eigenpy
4+
5+
dim = 100
6+
rng = np.random.default_rng()
7+
8+
A = rng.random((dim, dim))
9+
A = (A + A.T) * 0.5 + np.diag(5.0 + rng.random(dim))
10+
A = csc_matrix(A)
11+
12+
ilut = eigenpy.IncompleteLUT(A)
13+
assert ilut.info() == eigenpy.ComputationInfo.Success
14+
assert ilut.rows() == dim
15+
assert ilut.cols() == dim
16+
17+
X = rng.random((dim, 100))
18+
B = A.dot(X)
19+
X_est = ilut.solve(B)
20+
assert isinstance(X_est, np.ndarray)
21+
residual = np.linalg.norm(B - A.dot(X_est)) / np.linalg.norm(B)
22+
assert residual < 0.1
23+
24+
x = rng.random(dim)
25+
b = A.dot(x)
26+
x_est = ilut.solve(b)
27+
assert isinstance(x_est, np.ndarray)
28+
residual = np.linalg.norm(b - A.dot(x_est)) / np.linalg.norm(b)
29+
assert residual < 0.1
30+
31+
X_sparse = csc_matrix(rng.random((dim, 10)))
32+
B_sparse = A.dot(X_sparse).tocsc()
33+
if not B_sparse.has_sorted_indices:
34+
B_sparse.sort_indices()
35+
X_est_sparse = ilut.solve(B_sparse)
36+
assert isinstance(X_est_sparse, csc_matrix)
37+
38+
ilut.analyzePattern(A)
39+
ilut.factorize(A)
40+
assert ilut.info() == eigenpy.ComputationInfo.Success
41+
42+
ilut_params = eigenpy.IncompleteLUT(A, 1e-4, 15)
43+
assert ilut_params.info() == eigenpy.ComputationInfo.Success
44+
45+
ilut_set = eigenpy.IncompleteLUT()
46+
ilut_set.setDroptol(1e-3)
47+
ilut_set.setFillfactor(20)
48+
ilut_set.compute(A)
49+
assert ilut_set.info() == eigenpy.ComputationInfo.Success

0 commit comments

Comments
 (0)