diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0a38a1c3137b5c4f7627f8bebe3aa91fc77ead37..920fd3918a95dad11f3534271515c1f72f4998e7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,11 +6,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
## [Unreleased]
+### Added
+- Added id() helper to retrieve unique object identifier in Python ([#477](https://github.com/stack-of-tasks/eigenpy/pull/477))
+- Expose QR solvers ([#478](https://github.com/stack-of-tasks/eigenpy/pull/478))
+
## [3.6.0] - 2024-06-05
### Added
- Added a deprecation call policy shortcut ([#466](https://github.com/stack-of-tasks/eigenpy/pull/466))
-- Added id() helper to retrieve unique object identifier in Python ([#477](https://github.com/stack-of-tasks/eigenpy/pull/477))
### Fixed
- Fix register_symbolic_link_to_registered_type() for multiple successive registrations ([#471](https://github.com/stack-of-tasks/eigenpy/pull/471))
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 94bf5d85acf9716dc37d2e407748fc7de0aa9998..fa9e1e7b1e1e8bdf2fa169aa4276f5865fe69589 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -210,6 +210,11 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
include/eigenpy/decompositions/PermutationMatrix.hpp
include/eigenpy/decompositions/LDLT.hpp
include/eigenpy/decompositions/LLT.hpp
+ include/eigenpy/decompositions/QR.hpp
+ include/eigenpy/decompositions/HouseholderQR.hpp
+ include/eigenpy/decompositions/ColPivHouseholderQR.hpp
+ include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp
+ include/eigenpy/decompositions/FullPivHouseholderQR.hpp
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
include/eigenpy/decompositions/minres.hpp)
@@ -282,6 +287,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
src/decompositions/llt-solver.cpp
src/decompositions/ldlt-solver.cpp
src/decompositions/minres-solver.cpp
+ src/decompositions/qr-solvers.cpp
src/decompositions/eigen-solver.cpp
src/decompositions/self-adjoint-eigen-solver.cpp
src/decompositions/permutation-matrix.cpp
diff --git a/include/eigenpy/decompositions/ColPivHouseholderQR.hpp b/include/eigenpy/decompositions/ColPivHouseholderQR.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..95e1bd794e7f7122a6ffe6754d90ba0151c65ae6
--- /dev/null
+++ b/include/eigenpy/decompositions/ColPivHouseholderQR.hpp
@@ -0,0 +1,189 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_col_piv_houselholder_qr_hpp__
+#define __eigenpy_decompositions_col_piv_houselholder_qr_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+#include <Eigen/QR>
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct ColPivHouseholderQRSolverVisitor
+ : public boost::python::def_visitor<
+ ColPivHouseholderQRSolverVisitor<_MatrixType> > {
+ typedef _MatrixType MatrixType;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+ VectorXs;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+ MatrixType::Options>
+ MatrixXs;
+ typedef Eigen::ColPivHouseholderQR<MatrixType> Solver;
+ typedef Solver Self;
+
+ template <class PyClass>
+ void visit(PyClass &cl) const {
+ cl.def(bp::init<>(bp::arg("self"),
+ "Default constructor.\n"
+ "The default constructor is useful in cases in which the "
+ "user intends to perform decompositions via "
+ "HouseholderQR.compute(matrix)"))
+ .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
+ bp::args("self", "rows", "cols"),
+ "Default constructor with memory preallocation.\n"
+ "Like the default constructor but with preallocation of the "
+ "internal data according to the specified problem size. "))
+ .def(bp::init<MatrixType>(
+ bp::args("self", "matrix"),
+ "Constructs a QR factorization from a given matrix.\n"
+ "This constructor computes the QR factorization of the matrix "
+ "matrix by calling the method compute()."))
+
+ .def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
+ "Returns the absolute value of the determinant of the matrix of "
+ "which *this is the QR decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the QR decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices.")
+ .def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
+ "Returns the natural log of the absolute value of the determinant "
+ "of the matrix of which *this is the QR decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the QR decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices. This method is useful to "
+ "work around the risk of overflow/underflow that's inherent to "
+ "determinant computation.")
+ .def("dimensionOfKernel", &Self::dimensionOfKernel, bp::arg("self"),
+ "Returns the dimension of the kernel of the matrix of which *this "
+ "is the QR decomposition.")
+ .def("info", &Self::info, bp::arg("self"),
+ "Reports whether the QR factorization was successful.\n"
+ "Note: This function always returns Success. It is provided for "
+ "compatibility with other factorization routines.")
+ .def("isInjective", &Self::isInjective, bp::arg("self"),
+ "Returns true if the matrix associated with this QR decomposition "
+ "represents an injective linear map, i.e. has trivial kernel; "
+ "false otherwise.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("isInvertible", &Self::isInvertible, bp::arg("self"),
+ "Returns true if the matrix associated with the QR decomposition "
+ "is invertible.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("isSurjective", &Self::isSurjective, bp::arg("self"),
+ "Returns true if the matrix associated with this QR decomposition "
+ "represents a surjective linear map; false otherwise.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("maxPivot", &Self::maxPivot, bp::arg("self"),
+ "Returns the absolute value of the biggest pivot, i.e. the "
+ "biggest diagonal coefficient of U.")
+ .def("nonzeroPivots", &Self::nonzeroPivots, bp::arg("self"),
+ "Returns the number of nonzero pivots in the QR decomposition. "
+ "Here nonzero is meant in the exact sense, not in a fuzzy sense. "
+ "So that notion isn't really intrinsically interesting, but it is "
+ "still useful when implementing algorithms.")
+ .def("rank", &Self::rank, bp::arg("self"),
+ "Returns the rank of the matrix associated with the QR "
+ "decomposition.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+
+ .def("setThreshold",
+ (Self & (Self::*)(const RealScalar &)) & Self::setThreshold,
+ bp::args("self", "threshold"),
+ "Allows to prescribe a threshold to be used by certain methods, "
+ "such as rank(), who need to determine when pivots are to be "
+ "considered nonzero. This is not used for the QR decomposition "
+ "itself.\n"
+ "\n"
+ "When it needs to get the threshold value, Eigen calls "
+ "threshold(). By default, this uses a formula to automatically "
+ "determine a reasonable threshold. Once you have called the "
+ "present method setThreshold(const RealScalar&), your value is "
+ "used instead.\n"
+ "\n"
+ "Note: A pivot will be considered nonzero if its absolute value "
+ "is strictly greater than |pivot| ⩽ threshold×|maxpivot| where "
+ "maxpivot is the biggest pivot.",
+ bp::return_self<>())
+ .def("threshold", &Self::threshold, bp::arg("self"),
+ "Returns the threshold that will be used by certain methods such "
+ "as rank().")
+
+ .def("matrixQR", &Self::matrixQR, bp::arg("self"),
+ "Returns the matrix where the Householder QR decomposition is "
+ "stored in a LAPACK-compatible way.",
+ bp::return_value_policy<bp::copy_const_reference>())
+ .def("matrixR", &Self::matrixR, bp::arg("self"),
+ "Returns the matrix where the result Householder QR is stored.",
+ bp::return_value_policy<bp::copy_const_reference>())
+
+ .def(
+ "compute",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
+ Solver::compute,
+ bp::args("self", "matrix"),
+ "Computes the QR factorization of given matrix.",
+ bp::return_self<>())
+
+ .def("inverse", inverse, bp::arg("self"),
+ "Returns the inverse of the matrix associated with the QR "
+ "decomposition..")
+
+ .def("solve", &solve<MatrixXs>, 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.");
+ }
+
+ static void expose() {
+ static const std::string classname =
+ "ColPivHouseholderQR" + scalar_name<Scalar>::shortname();
+ expose(classname);
+ }
+
+ static void expose(const std::string &name) {
+ bp::class_<Solver>(
+ name.c_str(),
+ "This class performs a rank-revealing QR decomposition of a matrix A "
+ "into matrices P, Q and R such that:\n"
+ "AP=QR\n"
+ "by using Householder transformations. Here, P is a permutation "
+ "matrix, Q a unitary matrix and R an upper triangular matrix.\n"
+ "\n"
+ "This decomposition performs column pivoting in order to be "
+ "rank-revealing and improve numerical stability. It is slower than "
+ "HouseholderQR, and faster than FullPivHouseholderQR.",
+ bp::no_init)
+ .def(ColPivHouseholderQRSolverVisitor())
+ .def(IdVisitor<Solver>());
+ }
+
+ private:
+ template <typename MatrixOrVector>
+ static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+ return self.solve(vec);
+ }
+ static MatrixXs inverse(const Self &self) { return self.inverse(); }
+};
+
+} // namespace eigenpy
+
+#endif // ifndef __eigenpy_decompositions_col_piv_houselholder_qr_hpp__
diff --git a/include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp b/include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5cfc4a1c9a0399563128aa5c9f3c375d39be08a5
--- /dev/null
+++ b/include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp
@@ -0,0 +1,204 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_complete_orthogonal_decomposition_hpp__
+#define __eigenpy_decompositions_complete_orthogonal_decomposition_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+#include <Eigen/QR>
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct CompleteOrthogonalDecompositionSolverVisitor
+ : public boost::python::def_visitor<
+ CompleteOrthogonalDecompositionSolverVisitor<_MatrixType> > {
+ typedef _MatrixType MatrixType;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+ VectorXs;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+ MatrixType::Options>
+ MatrixXs;
+ typedef Eigen::CompleteOrthogonalDecomposition<MatrixType> Solver;
+ typedef Solver Self;
+
+ template <class PyClass>
+ void visit(PyClass &cl) const {
+ cl.def(bp::init<>(bp::arg("self"),
+ "Default constructor.\n"
+ "The default constructor is useful in cases in which the "
+ "user intends to perform decompositions via "
+ "CompleteOrthogonalDecomposition.compute(matrix)"))
+ .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
+ bp::args("self", "rows", "cols"),
+ "Default constructor with memory preallocation.\n"
+ "Like the default constructor but with preallocation of the "
+ "internal data according to the specified problem size. "))
+ .def(bp::init<MatrixType>(bp::args("self", "matrix"),
+ "Constructs a complete orthogonal "
+ "factorization from a given matrix.\n"
+ "This constructor computes the complete "
+ "orthogonal factorization of the matrix "
+ "matrix by calling the method compute()."))
+
+ .def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
+ "Returns the absolute value of the determinant of the matrix "
+ "associated with the complete orthogonal decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the complete orthogonal "
+ "decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices.")
+ .def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
+ "Returns the natural log of the absolute value of the determinant "
+ "of the matrix of which *this is the complete orthogonal "
+ "decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the complete orthogonal "
+ "decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices. This method is useful to "
+ "work around the risk of overflow/underflow that's inherent to "
+ "determinant computation.")
+ .def("dimensionOfKernel", &Self::dimensionOfKernel, bp::arg("self"),
+ "Returns the dimension of the kernel of the matrix of which *this "
+ "is the complete orthogonal decomposition.")
+ .def("info", &Self::info, bp::arg("self"),
+ "Reports whether the complete orthogonal factorization was "
+ "successful.\n"
+ "Note: This function always returns Success. It is provided for "
+ "compatibility with other factorization routines.")
+ .def("isInjective", &Self::isInjective, bp::arg("self"),
+ "Returns true if the matrix associated with this complete "
+ "orthogonal decomposition "
+ "represents an injective linear map, i.e. has trivial kernel; "
+ "false otherwise.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("isInvertible", &Self::isInvertible, bp::arg("self"),
+ "Returns true if the matrix associated with the complete "
+ "orthogonal decomposition "
+ "is invertible.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("isSurjective", &Self::isSurjective, bp::arg("self"),
+ "Returns true if the matrix associated with this complete "
+ "orthogonal decomposition "
+ "represents a surjective linear map; false otherwise.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("maxPivot", &Self::maxPivot, bp::arg("self"),
+ "Returns the absolute value of the biggest pivot, i.e. the "
+ "biggest diagonal coefficient of U.")
+ .def("nonzeroPivots", &Self::nonzeroPivots, bp::arg("self"),
+ "Returns the number of nonzero pivots in the complete orthogonal "
+ "decomposition. "
+ "Here nonzero is meant in the exact sense, not in a fuzzy sense. "
+ "So that notion isn't really intrinsically interesting, but it is "
+ "still useful when implementing algorithms.")
+ .def("rank", &Self::rank, bp::arg("self"),
+ "Returns the rank of the matrix associated with the complete "
+ "orthogonal "
+ "decomposition.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+
+ .def("setThreshold",
+ (Self & (Self::*)(const RealScalar &)) & Self::setThreshold,
+ bp::args("self", "threshold"),
+ "Allows to prescribe a threshold to be used by certain methods, "
+ "such as rank(), who need to determine when pivots are to be "
+ "considered nonzero. This is not used for the complete orthogonal "
+ "decomposition "
+ "itself.\n"
+ "\n"
+ "When it needs to get the threshold value, Eigen calls "
+ "threshold(). By default, this uses a formula to automatically "
+ "determine a reasonable threshold. Once you have called the "
+ "present method setThreshold(const RealScalar&), your value is "
+ "used instead.\n"
+ "\n"
+ "Note: A pivot will be considered nonzero if its absolute value "
+ "is strictly greater than |pivot| ⩽ threshold×|maxpivot| where "
+ "maxpivot is the biggest pivot.",
+ bp::return_self<>())
+ .def("threshold", &Self::threshold, bp::arg("self"),
+ "Returns the threshold that will be used by certain methods such "
+ "as rank().")
+
+ .def("matrixQTZ", &Self::matrixQTZ, bp::arg("self"),
+ "Returns the matrix where the complete orthogonal decomposition "
+ "is stored.",
+ bp::return_value_policy<bp::copy_const_reference>())
+ .def("matrixT", &Self::matrixT, bp::arg("self"),
+ "Returns the matrix where the complete orthogonal decomposition "
+ "is stored.",
+ bp::return_value_policy<bp::copy_const_reference>())
+ .def("matrixZ", &Self::matrixZ, bp::arg("self"),
+ "Returns the matrix Z.")
+
+ .def(
+ "compute",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
+ Solver::compute,
+ bp::args("self", "matrix"),
+ "Computes the complete orthogonal factorization of given matrix.",
+ bp::return_self<>())
+
+ .def("pseudoInverse", pseudoInverse, bp::arg("self"),
+ "Returns the pseudo-inverse of the matrix associated with the "
+ "complete orthogonal "
+ "decomposition.")
+
+ .def("solve", &solve<MatrixXs>, 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.");
+ }
+
+ static void expose() {
+ static const std::string classname =
+ "CompleteOrthogonalDecomposition" + scalar_name<Scalar>::shortname();
+ expose(classname);
+ }
+
+ static void expose(const std::string &name) {
+ bp::class_<Solver>(
+ name.c_str(),
+ "This class performs a rank-revealing complete orthogonal "
+ "decomposition of a matrix A into matrices P, Q, T, and Z such that:\n"
+ "AP=Q[T000]Z"
+ "by using Householder transformations. Here, P is a permutation "
+ "matrix, Q and Z are unitary matrices and T an upper triangular matrix "
+ "of size rank-by-rank. A may be rank deficient.",
+ bp::no_init)
+ .def(CompleteOrthogonalDecompositionSolverVisitor())
+ .def(IdVisitor<Solver>());
+ }
+
+ private:
+ template <typename MatrixOrVector>
+ static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+ return self.solve(vec);
+ }
+ static MatrixXs pseudoInverse(const Self &self) {
+ return self.pseudoInverse();
+ }
+};
+
+} // namespace eigenpy
+
+#endif // ifndef
+ // __eigenpy_decompositions_complete_orthogonal_decomposition_hpp__
diff --git a/include/eigenpy/decompositions/EigenSolver.hpp b/include/eigenpy/decompositions/EigenSolver.hpp
index 14583b2e3e9cade3d0bd2d42f7a9b8a162eaa824..5aa91cb8d2c05ebb10c1e483ceae7b3c40df6334 100644
--- a/include/eigenpy/decompositions/EigenSolver.hpp
+++ b/include/eigenpy/decompositions/EigenSolver.hpp
@@ -2,8 +2,8 @@
* Copyright 2020 INRIA
*/
-#ifndef __eigenpy_decomposition_eigen_solver_hpp__
-#define __eigenpy_decomposition_eigen_solver_hpp__
+#ifndef __eigenpy_decompositions_eigen_solver_hpp__
+#define __eigenpy_decompositions_eigen_solver_hpp__
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
@@ -90,4 +90,4 @@ struct EigenSolverVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_eigen_solver_hpp__
+#endif // ifndef __eigenpy_decompositions_eigen_solver_hpp__
diff --git a/include/eigenpy/decompositions/FullPivHouseholderQR.hpp b/include/eigenpy/decompositions/FullPivHouseholderQR.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..70dad1e454798167dc7b58bb5c75362294bd7d50
--- /dev/null
+++ b/include/eigenpy/decompositions/FullPivHouseholderQR.hpp
@@ -0,0 +1,183 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_full_piv_houselholder_qr_hpp__
+#define __eigenpy_decompositions_full_piv_houselholder_qr_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+#include <Eigen/QR>
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct FullPivHouseholderQRSolverVisitor
+ : public boost::python::def_visitor<
+ FullPivHouseholderQRSolverVisitor<_MatrixType> > {
+ typedef _MatrixType MatrixType;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+ VectorXs;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+ MatrixType::Options>
+ MatrixXs;
+ typedef Eigen::FullPivHouseholderQR<MatrixType> Solver;
+ typedef Solver Self;
+
+ template <class PyClass>
+ void visit(PyClass &cl) const {
+ cl.def(bp::init<>(bp::arg("self"),
+ "Default constructor.\n"
+ "The default constructor is useful in cases in which the "
+ "user intends to perform decompositions via "
+ "HouseholderQR.compute(matrix)"))
+ .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
+ bp::args("self", "rows", "cols"),
+ "Default constructor with memory preallocation.\n"
+ "Like the default constructor but with preallocation of the "
+ "internal data according to the specified problem size. "))
+ .def(bp::init<MatrixType>(
+ bp::args("self", "matrix"),
+ "Constructs a QR factorization from a given matrix.\n"
+ "This constructor computes the QR factorization of the matrix "
+ "matrix by calling the method compute()."))
+
+ .def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
+ "Returns the absolute value of the determinant of the matrix of "
+ "which *this is the QR decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the QR decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices.")
+ .def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
+ "Returns the natural log of the absolute value of the determinant "
+ "of the matrix of which *this is the QR decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the QR decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices. This method is useful to "
+ "work around the risk of overflow/underflow that's inherent to "
+ "determinant computation.")
+ .def("dimensionOfKernel", &Self::dimensionOfKernel, bp::arg("self"),
+ "Returns the dimension of the kernel of the matrix of which *this "
+ "is the QR decomposition.")
+ .def("isInjective", &Self::isInjective, bp::arg("self"),
+ "Returns true if the matrix associated with this QR decomposition "
+ "represents an injective linear map, i.e. has trivial kernel; "
+ "false otherwise.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("isInvertible", &Self::isInvertible, bp::arg("self"),
+ "Returns true if the matrix associated with the QR decomposition "
+ "is invertible.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("isSurjective", &Self::isSurjective, bp::arg("self"),
+ "Returns true if the matrix associated with this QR decomposition "
+ "represents a surjective linear map; false otherwise.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+ .def("maxPivot", &Self::maxPivot, bp::arg("self"),
+ "Returns the absolute value of the biggest pivot, i.e. the "
+ "biggest diagonal coefficient of U.")
+ .def("nonzeroPivots", &Self::nonzeroPivots, bp::arg("self"),
+ "Returns the number of nonzero pivots in the QR decomposition. "
+ "Here nonzero is meant in the exact sense, not in a fuzzy sense. "
+ "So that notion isn't really intrinsically interesting, but it is "
+ "still useful when implementing algorithms.")
+ .def("rank", &Self::rank, bp::arg("self"),
+ "Returns the rank of the matrix associated with the QR "
+ "decomposition.\n"
+ "\n"
+ "Note: This method has to determine which pivots should be "
+ "considered nonzero. For that, it uses the threshold value that "
+ "you can control by calling setThreshold(threshold).")
+
+ .def("setThreshold",
+ (Self & (Self::*)(const RealScalar &)) & Self::setThreshold,
+ bp::args("self", "threshold"),
+ "Allows to prescribe a threshold to be used by certain methods, "
+ "such as rank(), who need to determine when pivots are to be "
+ "considered nonzero. This is not used for the QR decomposition "
+ "itself.\n"
+ "\n"
+ "When it needs to get the threshold value, Eigen calls "
+ "threshold(). By default, this uses a formula to automatically "
+ "determine a reasonable threshold. Once you have called the "
+ "present method setThreshold(const RealScalar&), your value is "
+ "used instead.\n"
+ "\n"
+ "Note: A pivot will be considered nonzero if its absolute value "
+ "is strictly greater than |pivot| ⩽ threshold×|maxpivot| where "
+ "maxpivot is the biggest pivot.",
+ bp::return_self<>())
+ .def("threshold", &Self::threshold, bp::arg("self"),
+ "Returns the threshold that will be used by certain methods such "
+ "as rank().")
+
+ .def("matrixQR", &Self::matrixQR, bp::arg("self"),
+ "Returns the matrix where the Householder QR decomposition is "
+ "stored in a LAPACK-compatible way.",
+ bp::return_value_policy<bp::copy_const_reference>())
+
+ .def(
+ "compute",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
+ Solver::compute,
+ bp::args("self", "matrix"),
+ "Computes the QR factorization of given matrix.",
+ bp::return_self<>())
+
+ .def("inverse", inverse, bp::arg("self"),
+ "Returns the inverse of the matrix associated with the QR "
+ "decomposition..")
+
+ .def("solve", &solve<MatrixXs>, 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.");
+ }
+
+ static void expose() {
+ static const std::string classname =
+ "FullPivHouseholderQR" + scalar_name<Scalar>::shortname();
+ expose(classname);
+ }
+
+ static void expose(const std::string &name) {
+ bp::class_<Solver>(
+ name.c_str(),
+ "This class performs a rank-revealing QR decomposition of a matrix A "
+ "into matrices P, P', Q and R such that:\n"
+ "PAP′=QR\n"
+ "by using Householder transformations. Here, P and P' are permutation "
+ "matrices, Q a unitary matrix and R an upper triangular matrix.\n"
+ "\n"
+ "This decomposition performs a very prudent full pivoting in order to "
+ "be rank-revealing and achieve optimal numerical stability. The "
+ "trade-off is that it is slower than HouseholderQR and "
+ "ColPivHouseholderQR.",
+ bp::no_init)
+ .def(FullPivHouseholderQRSolverVisitor())
+ .def(IdVisitor<Solver>());
+ }
+
+ private:
+ template <typename MatrixOrVector>
+ static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+ return self.solve(vec);
+ }
+ static MatrixXs inverse(const Self &self) { return self.inverse(); }
+};
+
+} // namespace eigenpy
+
+#endif // ifndef __eigenpy_decompositions_full_piv_houselholder_qr_hpp__
diff --git a/include/eigenpy/decompositions/HouseholderQR.hpp b/include/eigenpy/decompositions/HouseholderQR.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..895000f4a0ac7131762f5ec0e42f8dbd23845ee5
--- /dev/null
+++ b/include/eigenpy/decompositions/HouseholderQR.hpp
@@ -0,0 +1,118 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_houselholder_qr_hpp__
+#define __eigenpy_decompositions_houselholder_qr_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+
+#include <Eigen/QR>
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct HouseholderQRSolverVisitor
+ : public boost::python::def_visitor<
+ HouseholderQRSolverVisitor<_MatrixType> > {
+ typedef _MatrixType MatrixType;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+ VectorXs;
+ typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+ MatrixType::Options>
+ MatrixXs;
+ typedef Eigen::HouseholderQR<MatrixType> Solver;
+ typedef Solver Self;
+
+ template <class PyClass>
+ void visit(PyClass &cl) const {
+ cl.def(bp::init<>(bp::arg("self"),
+ "Default constructor.\n"
+ "The default constructor is useful in cases in which the "
+ "user intends to perform decompositions via "
+ "HouseholderQR.compute(matrix)"))
+ .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
+ bp::args("self", "rows", "cols"),
+ "Default constructor with memory preallocation.\n"
+ "Like the default constructor but with preallocation of the "
+ "internal data according to the specified problem size. "))
+ .def(bp::init<MatrixType>(
+ bp::args("self", "matrix"),
+ "Constructs a QR factorization from a given matrix.\n"
+ "This constructor computes the QR factorization of the matrix "
+ "matrix by calling the method compute()."))
+
+ .def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
+ "Returns the absolute value of the determinant of the matrix of "
+ "which *this is the QR decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the QR decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices.")
+ .def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
+ "Returns the natural log of the absolute value of the determinant "
+ "of the matrix of which *this is the QR decomposition.\n"
+ "It has only linear complexity (that is, O(n) where n is the "
+ "dimension of the square matrix) as the QR decomposition has "
+ "already been computed.\n"
+ "Note: This is only for square matrices. This method is useful to "
+ "work around the risk of overflow/underflow that's inherent to "
+ "determinant computation.")
+
+ .def("matrixQR", &Self::matrixQR, bp::arg("self"),
+ "Returns the matrix where the Householder QR decomposition is "
+ "stored in a LAPACK-compatible way.",
+ bp::return_value_policy<bp::copy_const_reference>())
+
+ .def(
+ "compute",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
+ Solver::compute,
+ bp::args("self", "matrix"),
+ "Computes the QR factorization of given matrix.",
+ bp::return_self<>())
+
+ .def("solve", &solve<MatrixXs>, 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.");
+ }
+
+ static void expose() {
+ static const std::string classname =
+ "HouseholderQR" + scalar_name<Scalar>::shortname();
+ expose(classname);
+ }
+
+ static void expose(const std::string &name) {
+ bp::class_<Solver>(
+ name.c_str(),
+ "This class performs a QR decomposition of a matrix A into matrices Q "
+ "and R such that A=QR by using Householder transformations.\n"
+ "Here, Q a unitary matrix and R an upper triangular matrix. The result "
+ "is stored in a compact way compatible with LAPACK.\n"
+ "\n"
+ "Note that no pivoting is performed. This is not a rank-revealing "
+ "decomposition. If you want that feature, use FullPivHouseholderQR or "
+ "ColPivHouseholderQR instead.\n"
+ "\n"
+ "This Householder QR decomposition is faster, but less numerically "
+ "stable and less feature-full than FullPivHouseholderQR or "
+ "ColPivHouseholderQR.",
+ bp::no_init)
+ .def(HouseholderQRSolverVisitor())
+ .def(IdVisitor<Solver>());
+ }
+
+ private:
+ template <typename MatrixOrVector>
+ static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+ return self.solve(vec);
+ }
+};
+
+} // namespace eigenpy
+
+#endif // ifndef __eigenpy_decompositions_houselholder_qr_hpp__
diff --git a/include/eigenpy/decompositions/LDLT.hpp b/include/eigenpy/decompositions/LDLT.hpp
index b148a23a2ddbc652cdb7a75e743c632b33a4931c..4f67d2713269760358a7088afca3a1708703a017 100644
--- a/include/eigenpy/decompositions/LDLT.hpp
+++ b/include/eigenpy/decompositions/LDLT.hpp
@@ -2,8 +2,8 @@
* Copyright 2020-2024 INRIA
*/
-#ifndef __eigenpy_decomposition_ldlt_hpp__
-#define __eigenpy_decomposition_ldlt_hpp__
+#ifndef __eigenpy_decompositions_ldlt_hpp__
+#define __eigenpy_decompositions_ldlt_hpp__
#include <Eigen/Cholesky>
#include <Eigen/Core>
@@ -141,4 +141,4 @@ struct LDLTSolverVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_ldlt_hpp__
+#endif // ifndef __eigenpy_decompositions_ldlt_hpp__
diff --git a/include/eigenpy/decompositions/LLT.hpp b/include/eigenpy/decompositions/LLT.hpp
index 695c37300cc2072fe2ca7f96a26c9ce8aecbf85a..abf898ec3f672e90e033698bce19ee27ea0be072 100644
--- a/include/eigenpy/decompositions/LLT.hpp
+++ b/include/eigenpy/decompositions/LLT.hpp
@@ -2,8 +2,8 @@
* Copyright 2020-2024 INRIA
*/
-#ifndef __eigenpy_decomposition_llt_hpp__
-#define __eigenpy_decomposition_llt_hpp__
+#ifndef __eigenpy_decompositions_llt_hpp__
+#define __eigenpy_decompositions_llt_hpp__
#include <Eigen/Cholesky>
#include <Eigen/Core>
@@ -131,4 +131,4 @@ struct LLTSolverVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_llt_hpp__
+#endif // ifndef __eigenpy_decompositions_llt_hpp__
diff --git a/include/eigenpy/decompositions/PermutationMatrix.hpp b/include/eigenpy/decompositions/PermutationMatrix.hpp
index 094bdbaf1144c0258c9e2ba36aa4910550f4ab72..d66396c85f8a0f1d796fa0320ae35887666bd207 100644
--- a/include/eigenpy/decompositions/PermutationMatrix.hpp
+++ b/include/eigenpy/decompositions/PermutationMatrix.hpp
@@ -2,8 +2,8 @@
* Copyright 2024 INRIA
*/
-#ifndef __eigenpy_decomposition_permutation_matrix_hpp__
-#define __eigenpy_decomposition_permutation_matrix_hpp__
+#ifndef __eigenpy_decompositions_permutation_matrix_hpp__
+#define __eigenpy_decompositions_permutation_matrix_hpp__
#include "eigenpy/eigenpy.hpp"
#include "eigenpy/eigen/EigenBase.hpp"
@@ -104,4 +104,4 @@ struct PermutationMatrixVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_permutation_matrix_hpp__
+#endif // ifndef __eigenpy_decompositions_permutation_matrix_hpp__
diff --git a/include/eigenpy/decompositions/QR.hpp b/include/eigenpy/decompositions/QR.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..14868a1fc262b86bb80cddd5de14ea8b35d3b066
--- /dev/null
+++ b/include/eigenpy/decompositions/QR.hpp
@@ -0,0 +1,13 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_qr_hpp__
+#define __eigenpy_decompositions_qr_hpp__
+
+#include "eigenpy/decompositions/HouseholderQR.hpp"
+#include "eigenpy/decompositions/FullPivHouseholderQR.hpp"
+#include "eigenpy/decompositions/ColPivHouseholderQR.hpp"
+#include "eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp"
+
+#endif // ifndef __eigenpy_decompositions_qr_hpp__
diff --git a/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp b/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
index b2587a802459c192b6358a6abf3d4daf75292828..8a378fa3c88592847e1ccae377fe171baf28c858 100644
--- a/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
+++ b/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
@@ -2,8 +2,8 @@
* Copyright 2020-2024 INRIA
*/
-#ifndef __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
-#define __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
+#ifndef __eigenpy_decompositions_self_adjoint_eigen_solver_hpp__
+#define __eigenpy_decompositions_self_adjoint_eigen_solver_hpp__
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
@@ -102,4 +102,4 @@ struct SelfAdjointEigenSolverVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
+#endif // ifndef __eigenpy_decompositions_self_adjoint_eigen_solver_hpp__
diff --git a/include/eigenpy/decompositions/minres.hpp b/include/eigenpy/decompositions/minres.hpp
index 7f5019c4b61f708f546c6cdbfc0638305fa6e65c..03ee7d19b260bba67608ca8f59c79dc403d58545 100644
--- a/include/eigenpy/decompositions/minres.hpp
+++ b/include/eigenpy/decompositions/minres.hpp
@@ -2,8 +2,8 @@
* Copyright 2021 INRIA
*/
-#ifndef __eigenpy_decomposition_minres_hpp__
-#define __eigenpy_decomposition_minres_hpp__
+#ifndef __eigenpy_decompositions_minres_hpp__
+#define __eigenpy_decompositions_minres_hpp__
#include <Eigen/Core>
#include <iostream>
@@ -178,4 +178,4 @@ struct MINRESSolverVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_minres_hpp__
+#endif // ifndef __eigenpy_decompositions_minres_hpp__
diff --git a/include/eigenpy/decompositions/sparse/LDLT.hpp b/include/eigenpy/decompositions/sparse/LDLT.hpp
index 0d552e4e86cfb2f7c2b342d6e0ddbd34a826dac6..53772aede9ed9b292ff235d0d3371957524c4f5f 100644
--- a/include/eigenpy/decompositions/sparse/LDLT.hpp
+++ b/include/eigenpy/decompositions/sparse/LDLT.hpp
@@ -2,8 +2,8 @@
* Copyright 2024 INRIA
*/
-#ifndef __eigenpy_decomposition_sparse_ldlt_hpp__
-#define __eigenpy_decomposition_sparse_ldlt_hpp__
+#ifndef __eigenpy_decompositions_sparse_ldlt_hpp__
+#define __eigenpy_decompositions_sparse_ldlt_hpp__
#include "eigenpy/eigenpy.hpp"
#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
@@ -69,4 +69,4 @@ struct SimplicialLDLTVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_sparse_ldlt_hpp__
+#endif // ifndef __eigenpy_decompositions_sparse_ldlt_hpp__
diff --git a/include/eigenpy/decompositions/sparse/LLT.hpp b/include/eigenpy/decompositions/sparse/LLT.hpp
index 3d3fa0d562e99fc468bfb976cd19dfe0e4fd8d4e..4e521636547f7a51d42bc40a5ba33fdee7ddc4fd 100644
--- a/include/eigenpy/decompositions/sparse/LLT.hpp
+++ b/include/eigenpy/decompositions/sparse/LLT.hpp
@@ -2,8 +2,8 @@
* Copyright 2024 INRIA
*/
-#ifndef __eigenpy_decomposition_sparse_llt_hpp__
-#define __eigenpy_decomposition_sparse_llt_hpp__
+#ifndef __eigenpy_decompositions_sparse_llt_hpp__
+#define __eigenpy_decompositions_sparse_llt_hpp__
#include "eigenpy/eigenpy.hpp"
#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
@@ -64,4 +64,4 @@ struct SimplicialLLTVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_sparse_llt_hpp__
+#endif // ifndef __eigenpy_decompositions_sparse_llt_hpp__
diff --git a/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp b/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
index ad109140b08eeb79e8705ce419086203bd6a6721..912674b5f02a0a8d57f44ce601f8996a3469aae1 100644
--- a/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
+++ b/include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
@@ -2,8 +2,8 @@
* Copyright 2024 INRIA
*/
-#ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
-#define __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+#ifndef __eigenpy_decompositions_sparse_simplicial_cholesky_hpp__
+#define __eigenpy_decompositions_sparse_simplicial_cholesky_hpp__
#include "eigenpy/eigenpy.hpp"
#include "eigenpy/eigen/EigenBase.hpp"
@@ -91,4 +91,4 @@ struct SimplicialCholeskyVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+#endif // ifndef __eigenpy_decompositions_sparse_simplicial_cholesky_hpp__
diff --git a/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp b/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
index 2afc1b917d1060ef037f5186f3b885f340cfc87d..a1b3c8522edc2fdefb318eb309d557f39488f4ba 100644
--- a/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
+++ b/include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
@@ -2,8 +2,8 @@
* Copyright 2024 INRIA
*/
-#ifndef __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
-#define __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
+#ifndef __eigenpy_decompositions_sparse_sparse_solver_base_hpp__
+#define __eigenpy_decompositions_sparse_sparse_solver_base_hpp__
#include "eigenpy/eigenpy.hpp"
#include "eigenpy/eigen/EigenBase.hpp"
@@ -51,4 +51,4 @@ struct SparseSolverBaseVisitor
} // namespace eigenpy
-#endif // ifndef __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
+#endif // ifndef __eigenpy_decompositions_sparse_sparse_solver_base_hpp__
diff --git a/src/decompositions/decompositions.cpp b/src/decompositions/decompositions.cpp
index a497db1e8012c5ba0efd6058a77c01942405289b..77ae22f295b6097420454c60f91de861ad7e2707 100644
--- a/src/decompositions/decompositions.cpp
+++ b/src/decompositions/decompositions.cpp
@@ -12,6 +12,7 @@ void exposeEigenSolver();
void exposeSelfAdjointEigenSolver();
void exposeLLTSolver();
void exposeLDLTSolver();
+void exposeQRSolvers();
void exposeMINRESSolver();
void exposeSimplicialLLTSolver();
void exposeSimplicialLDLTSolver();
@@ -24,6 +25,7 @@ void exposeDecompositions() {
exposeSelfAdjointEigenSolver();
exposeLLTSolver();
exposeLDLTSolver();
+ exposeQRSolvers();
exposeMINRESSolver();
{
diff --git a/src/decompositions/qr-solvers.cpp b/src/decompositions/qr-solvers.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..75ab771fa4c87f6a518b916346532fac42cee571
--- /dev/null
+++ b/src/decompositions/qr-solvers.cpp
@@ -0,0 +1,16 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#include "eigenpy/decompositions/QR.hpp"
+
+namespace eigenpy {
+void exposeQRSolvers() {
+ using namespace Eigen;
+ HouseholderQRSolverVisitor<MatrixXd>::expose("HouseholderQR");
+ FullPivHouseholderQRSolverVisitor<MatrixXd>::expose("FullPivHouseholderQR");
+ ColPivHouseholderQRSolverVisitor<MatrixXd>::expose("ColPivHouseholderQR");
+ CompleteOrthogonalDecompositionSolverVisitor<MatrixXd>::expose(
+ "CompleteOrthogonalDecomposition");
+}
+} // namespace eigenpy
diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt
index 38f75745f9f7672bf953bb8bba2ccd36190c6059..8d6571d1bade3aec7eced4e65c7feaeb02185771 100644
--- a/unittest/CMakeLists.txt
+++ b/unittest/CMakeLists.txt
@@ -136,6 +136,8 @@ add_python_eigenpy_lib_unit_test("py-LDLT" "unittest/python/test_LDLT.py")
add_python_eigenpy_lib_unit_test("py-id" "unittest/python/test_id.py")
+add_python_eigenpy_lib_unit_test("py-QR" "unittest/python/test_QR.py")
+
if(NOT WIN32)
add_python_eigenpy_lib_unit_test("py-MINRES" "unittest/python/test_MINRES.py")
endif(NOT WIN32)
diff --git a/unittest/python/test_QR.py b/unittest/python/test_QR.py
new file mode 100644
index 0000000000000000000000000000000000000000..aaf13cf09747f564fcc8d82c5b496de4e163d22e
--- /dev/null
+++ b/unittest/python/test_QR.py
@@ -0,0 +1,76 @@
+import numpy as np
+
+import eigenpy
+
+rows = 20
+cols = 100
+rng = np.random.default_rng()
+
+A = rng.random((rows, cols))
+
+# Test HouseholderQR decomposition
+householder_qr = eigenpy.HouseholderQR()
+householder_qr = eigenpy.HouseholderQR(rows, cols)
+householder_qr = eigenpy.HouseholderQR(A)
+
+householder_qr_eye = eigenpy.HouseholderQR(np.eye(rows, rows))
+X = rng.random((rows, 20))
+assert householder_qr_eye.absDeterminant() == 1.0
+assert householder_qr_eye.logAbsDeterminant() == 0.0
+
+Y = householder_qr_eye.solve(X)
+assert (X == Y).all()
+
+# Test FullPivHouseholderQR decomposition
+fullpiv_householder_qr = eigenpy.FullPivHouseholderQR()
+fullpiv_householder_qr = eigenpy.FullPivHouseholderQR(rows, cols)
+fullpiv_householder_qr = eigenpy.FullPivHouseholderQR(A)
+
+fullpiv_householder_qr = eigenpy.FullPivHouseholderQR(np.eye(rows, rows))
+X = rng.random((rows, 20))
+assert fullpiv_householder_qr.absDeterminant() == 1.0
+assert fullpiv_householder_qr.logAbsDeterminant() == 0.0
+
+Y = fullpiv_householder_qr.solve(X)
+assert (X == Y).all()
+assert fullpiv_householder_qr.rank() == rows
+
+fullpiv_householder_qr.setThreshold(1e-8)
+assert fullpiv_householder_qr.threshold() == 1e-8
+assert eigenpy.is_approx(np.eye(rows, rows), fullpiv_householder_qr.inverse())
+
+# Test ColPivHouseholderQR decomposition
+colpiv_householder_qr = eigenpy.ColPivHouseholderQR()
+colpiv_householder_qr = eigenpy.ColPivHouseholderQR(rows, cols)
+colpiv_householder_qr = eigenpy.ColPivHouseholderQR(A)
+
+colpiv_householder_qr = eigenpy.ColPivHouseholderQR(np.eye(rows, rows))
+X = rng.random((rows, 20))
+assert colpiv_householder_qr.absDeterminant() == 1.0
+assert colpiv_householder_qr.logAbsDeterminant() == 0.0
+
+Y = colpiv_householder_qr.solve(X)
+assert (X == Y).all()
+assert colpiv_householder_qr.rank() == rows
+
+colpiv_householder_qr.setThreshold(1e-8)
+assert colpiv_householder_qr.threshold() == 1e-8
+assert eigenpy.is_approx(np.eye(rows, rows), colpiv_householder_qr.inverse())
+
+# Test CompleteOrthogonalDecomposition
+cod = eigenpy.CompleteOrthogonalDecomposition()
+cod = eigenpy.CompleteOrthogonalDecomposition(rows, cols)
+cod = eigenpy.CompleteOrthogonalDecomposition(A)
+
+cod = eigenpy.CompleteOrthogonalDecomposition(np.eye(rows, rows))
+X = rng.random((rows, 20))
+assert cod.absDeterminant() == 1.0
+assert cod.logAbsDeterminant() == 0.0
+
+Y = cod.solve(X)
+assert (X == Y).all()
+assert cod.rank() == rows
+
+cod.setThreshold(1e-8)
+assert cod.threshold() == 1e-8
+assert eigenpy.is_approx(np.eye(rows, rows), cod.pseudoInverse())