Skip to content

Commit 5f02e30

Browse files
Lucas-HaubertLucas-Haubert
committed
eigenvalues: ComplexEigenSolver
1 parent f70b8d7 commit 5f02e30

File tree

7 files changed

+169
-33
lines changed

7 files changed

+169
-33
lines changed

‎CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -236,6 +236,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
236236
${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS}
237237
include/eigenpy/decompositions/decompositions.hpp
238238
include/eigenpy/decompositions/EigenSolver.hpp
239+
include/eigenpy/decompositions/ComplexEigenSolver.hpp
239240
include/eigenpy/decompositions/PermutationMatrix.hpp
240241
include/eigenpy/decompositions/LDLT.hpp
241242
include/eigenpy/decompositions/LLT.hpp
@@ -318,6 +319,7 @@ set(${PROJECT_NAME}_SOLVERS_SOURCES src/solvers/preconditioners.cpp
318319
set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
319320
src/decompositions/decompositions.cpp
320321
src/decompositions/eigen-solver.cpp
322+
src/decompositions/complex-eigen-solver.cpp
321323
src/decompositions/llt-solver.cpp
322324
src/decompositions/ldlt-solver.cpp
323325
src/decompositions/bdcsvd-solver.cpp
@@ -328,7 +330,6 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
328330
src/decompositions/sparse-lu-solver.cpp
329331
src/decompositions/sparse-qr-solver.cpp
330332
src/decompositions/qr-solvers.cpp
331-
src/decompositions/eigen-solver.cpp
332333
src/decompositions/self-adjoint-eigen-solver.cpp
333334
src/decompositions/permutation-matrix.cpp
334335
src/decompositions/simplicial-llt-solver.cpp

‎include/eigenpy/decompositions/ComplexEigenSolver.hpp‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
/*
2+
* Copyright 2020 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_complex_eigen_solver_hpp__
6+
#define __eigenpy_decompositions_complex_eigen_solver_hpp__
7+
8+
#include <Eigen/Core>
9+
#include <Eigen/Eigenvalues>
10+
11+
#include "eigenpy/eigen-to-python.hpp"
12+
#include "eigenpy/eigenpy.hpp"
13+
#include "eigenpy/utils/scalar-name.hpp"
14+
15+
namespace eigenpy {
16+
17+
template <typename _MatrixType>
18+
struct ComplexEigenSolverVisitor : public boost::python::def_visitor<
19+
ComplexEigenSolverVisitor<_MatrixType>> {
20+
typedef _MatrixType MatrixType;
21+
typedef typename MatrixType::Scalar Scalar;
22+
typedef Eigen::ComplexEigenSolver<MatrixType> Solver;
23+
24+
template <class PyClass>
25+
void visit(PyClass& cl) const {
26+
cl.def(bp::init<>("Default constructor"))
27+
.def(bp::init<Eigen::DenseIndex>(
28+
bp::arg("size"), "Default constructor with memory preallocation"))
29+
.def(bp::init<MatrixType, bp::optional<bool>>(
30+
bp::args("matrix", "compute_eigen_vectors"),
31+
"Computes eigendecomposition of given matrix"))
32+
33+
.def("eigenvalues", &Solver::eigenvalues, bp::arg("self"),
34+
"Returns the eigenvalues of given matrix.",
35+
bp::return_internal_reference<>())
36+
.def("eigenvectors", &Solver::eigenvectors, bp::arg("self"),
37+
"Returns the eigenvectors of given matrix.",
38+
bp::return_internal_reference<>())
39+
40+
.def("compute", &ComplexEigenSolverVisitor::compute_proxy<MatrixType>,
41+
bp::args("self", "matrix"),
42+
"Computes the eigendecomposition of given matrix.",
43+
bp::return_self<>())
44+
.def("compute",
45+
(Solver &
46+
(Solver::*)(const Eigen::EigenBase<MatrixType>& matrix, bool)) &
47+
Solver::compute,
48+
bp::args("self", "matrix", "compute_eigen_vectors"),
49+
"Computes the eigendecomposition of given matrix.",
50+
bp::return_self<>())
51+
52+
.def("info", &Solver::info, bp::arg("self"),
53+
"NumericalIssue if the input contains INF or NaN values or "
54+
"overflow occured. Returns Success otherwise.")
55+
56+
.def("getMaxIterations", &Solver::getMaxIterations, bp::arg("self"),
57+
"Returns the maximum number of iterations.")
58+
.def("setMaxIterations", &Solver::setMaxIterations,
59+
bp::args("self", "max_iter"),
60+
"Sets the maximum number of iterations allowed.",
61+
bp::return_self<>());
62+
}
63+
64+
static void expose() {
65+
static const std::string classname =
66+
"ComplexEigenSolver" + scalar_name<Scalar>::shortname();
67+
expose(classname);
68+
}
69+
70+
static void expose(const std::string& name) {
71+
bp::class_<Solver>(name.c_str(), bp::no_init)
72+
.def(ComplexEigenSolverVisitor())
73+
.def(IdVisitor<Solver>());
74+
}
75+
76+
private:
77+
template <typename MatrixType>
78+
static Solver& compute_proxy(Solver& self,
79+
const Eigen::EigenBase<MatrixType>& matrix) {
80+
return self.compute(matrix);
81+
}
82+
};
83+
84+
} // namespace eigenpy
85+
86+
#endif // ifndef __eigenpy_decompositions_complex_eigen_solver_hpp__

‎include/eigenpy/decompositions/JacobiSVD.hpp‎

Lines changed: 48 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,6 @@ namespace eigenpy {
1818
template <typename _MatrixType>
1919
struct JacobiSVDVisitor
2020
: public boost::python::def_visitor<JacobiSVDVisitor<_MatrixType>> {
21-
2221
typedef _MatrixType MatrixType;
2322
typedef Eigen::JacobiSVD<MatrixType> Solver;
2423
typedef typename MatrixType::Scalar Scalar;
@@ -32,34 +31,39 @@ struct JacobiSVDVisitor
3231
.def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex, unsigned int>(
3332
bp::args("self", "rows", "cols", "computationOptions "),
3433
"Default Constructor with memory preallocation. \n\n"
35-
"Like the default constructor but with preallocation of the internal "
36-
"data according to the specified problem size and the computationOptions. "))
34+
"Like the default constructor but with preallocation of the "
35+
"internal "
36+
"data according to the specified problem size and the "
37+
"computationOptions. "))
3738
.def(bp::init<MatrixType>(
3839
bp::args("self", "matrix"),
3940
"Constructor performing the decomposition of given matrix. "))
4041
.def(bp::init<MatrixType, unsigned int>(
4142
bp::args("self", "matrix", "computationOptions "),
4243
"Constructor performing the decomposition of given matrix. \n\n"
43-
"One cannot request unitiaries using both the Options template parameter "
44-
"and the constructor. If possible, prefer using the Options template parameter."))
44+
"One cannot request unitiaries using both the Options template "
45+
"parameter "
46+
"and the constructor. If possible, prefer using the Options "
47+
"template parameter."))
4548

4649
.def("cols", &Solver::cols, bp::arg("self"),
4750
"Returns the number of columns. ")
48-
.def(
49-
"compute",
50-
(Solver & (Solver::*)(const MatrixType &matrix)) &
51-
Solver::compute,
52-
bp::args("self", "matrix"),
53-
"Method performing the decomposition of given matrix. Computes Thin/Full "
54-
"unitaries U/V if specified using the Options template parameter or the class constructor. ",
55-
bp::return_self<>())
56-
.def(
57-
"compute",
58-
(Solver & (Solver::*)(const MatrixType &matrix, unsigned int computationOptions)) &
59-
Solver::compute,
60-
bp::args("self", "matrix"),
61-
"Method performing the decomposition of given matrix, as specified by the computationOptions parameter. ",
62-
bp::return_self<>())
51+
.def("compute",
52+
(Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute,
53+
bp::args("self", "matrix"),
54+
"Method performing the decomposition of given matrix. Computes "
55+
"Thin/Full "
56+
"unitaries U/V if specified using the Options template parameter "
57+
"or the class constructor. ",
58+
bp::return_self<>())
59+
.def("compute",
60+
(Solver & (Solver::*)(const MatrixType &matrix,
61+
unsigned int computationOptions)) &
62+
Solver::compute,
63+
bp::args("self", "matrix"),
64+
"Method performing the decomposition of given matrix, as "
65+
"specified by the computationOptions parameter. ",
66+
bp::return_self<>())
6367
.def("rows", &Solver::rows, bp::arg("self"),
6468
"Returns the number of rows. . ")
6569

@@ -76,20 +80,32 @@ struct JacobiSVDVisitor
7680
bp::class_<Solver, boost::noncopyable>(
7781
name.c_str(),
7882
"Two-sided Jacobi SVD decomposition of a rectangular matrix. \n\n"
79-
"SVD decomposition consists in decomposing any n-by-p matrix A as a product "
80-
"A=USV∗ where U is a n-by-n unitary, V is a p-by-p unitary, and S is a n-by-p r"
81-
"eal positive matrix which is zero outside of its main diagonal; the diagonal "
82-
"entries of S are known as the singular values of A and the columns of U and V "
83-
"are known as the left and right singular vectors of A respectively. \n\n"
83+
"SVD decomposition consists in decomposing any n-by-p matrix A as a "
84+
"product "
85+
"A=USV∗ where U is a n-by-n unitary, V is a p-by-p unitary, and S is a "
86+
"n-by-p r"
87+
"eal positive matrix which is zero outside of its main diagonal; the "
88+
"diagonal "
89+
"entries of S are known as the singular values of A and the columns of "
90+
"U and V "
91+
"are known as the left and right singular vectors of A respectively. "
92+
"\n\n"
8493
"Singular values are always sorted in decreasing order. \n\n "
85-
"This JacobiSVD decomposition computes only the singular values by default. "
94+
"This JacobiSVD decomposition computes only the singular values by "
95+
"default. "
8696
"If you want U or V, you need to ask for them explicitly. \n\n"
87-
"You can ask for only thin U or V to be computed, meaning the following. "
88-
"In case of a rectangular n-by-p matrix, letting m be the smaller value among "
89-
"n and p, there are only m singular vectors; the remaining columns of U and V "
90-
"do not correspond to actual singular vectors. Asking for thin U or V means asking "
91-
"for only their m first columns to be formed. So U is then a n-by-m matrix, and V "
92-
"is then a p-by-m matrix. Notice that thin U and V are all you need for (least squares) "
97+
"You can ask for only thin U or V to be computed, meaning the "
98+
"following. "
99+
"In case of a rectangular n-by-p matrix, letting m be the smaller "
100+
"value among "
101+
"n and p, there are only m singular vectors; the remaining columns of "
102+
"U and V "
103+
"do not correspond to actual singular vectors. Asking for thin U or V "
104+
"means asking "
105+
"for only their m first columns to be formed. So U is then a n-by-m "
106+
"matrix, and V "
107+
"is then a p-by-m matrix. Notice that thin U and V are all you need "
108+
"for (least squares) "
93109
"solving.",
94110
bp::no_init)
95111
.def(JacobiSVDVisitor())

‎src/decompositions/complex-eigen-solver.cpp‎

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

‎src/decompositions/decompositions.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
namespace eigenpy {
1010

1111
void exposeEigenSolver();
12+
void exposeComplexEigenSolver();
1213
void exposeSelfAdjointEigenSolver();
1314
void exposeLLTSolver();
1415
void exposeLDLTSolver();
@@ -28,6 +29,7 @@ void exposeDecompositions() {
2829
using namespace Eigen;
2930

3031
exposeEigenSolver();
32+
exposeComplexEigenSolver();
3133
exposeSelfAdjointEigenSolver();
3234
exposeLLTSolver();
3335
exposeLDLTSolver();

‎unittest/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,9 @@ add_python_eigenpy_lib_unit_test("py-version" "unittest/python/test_version.py")
137137
add_python_eigenpy_lib_unit_test("py-eigen-solver"
138138
"unittest/python/test_eigen_solver.py")
139139

140+
add_python_eigenpy_lib_unit_test("py-complex-eigen-solver"
141+
"unittest/python/test_complex_eigen_solver.py")
142+
140143
add_python_eigenpy_lib_unit_test(
141144
"py-self-adjoint-eigen-solver"
142145
"unittest/python/test_self_adjoint_eigen_solver.py")

‎unittest/python/test_complex_eigen_solver.py‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
import numpy as np
2+
3+
import eigenpy
4+
5+
dim = 100
6+
rng = np.random.default_rng()
7+
A = rng.random((dim, dim))
8+
9+
es = eigenpy.ComplexEigenSolver(A)
10+
11+
V = es.eigenvectors()
12+
D = es.eigenvalues()
13+
14+
assert eigenpy.is_approx(A.dot(V).real, V.dot(np.diag(D)).real)
15+
assert eigenpy.is_approx(A.dot(V).imag, V.dot(np.diag(D)).imag)

0 commit comments

Comments
 (0)