diff --git a/CHANGELOG.md b/CHANGELOG.md index bd14ced6f..a48e365d9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). ### Added +- Add additional decompositions and solvers from Eigen ([#571](https://github.com/stack-of-tasks/eigenpy/pull/571)) + +### Added + - Docker images `ghcr.io/stack-of-tasks/eigenpy` ([#575](https://github.com/stack-of-tasks/eigenpy/pull/575)) ### Changed diff --git a/CMakeLists.txt b/CMakeLists.txt index e82a3c6ae..e767907c8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -190,10 +190,14 @@ set(${PROJECT_NAME}_SOLVERS_HEADERS include/eigenpy/solvers/preconditioners.hpp include/eigenpy/solvers/IterativeSolverBase.hpp include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp + include/eigenpy/solvers/BiCGSTAB.hpp + include/eigenpy/solvers/MINRES.hpp include/eigenpy/solvers/ConjugateGradient.hpp include/eigenpy/solvers/SparseSolverBase.hpp include/eigenpy/solvers/BasicPreconditioners.hpp - include/eigenpy/solvers/BFGSPreconditioners.hpp) + include/eigenpy/solvers/BFGSPreconditioners.hpp + include/eigenpy/solvers/IncompleteCholesky.hpp + include/eigenpy/solvers/IncompleteLUT.hpp) set(${PROJECT_NAME}_EIGEN_HEADERS include/eigenpy/eigen/EigenBase.hpp) @@ -205,13 +209,18 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_CHOLMOD_HEADERS include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp) set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_ACCELERATE_HEADERS + include/eigenpy/decompositions/sparse/accelerate/Accelerate.hpp include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp) set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS - include/eigenpy/decompositions/sparse/LLT.hpp - include/eigenpy/decompositions/sparse/LDLT.hpp + include/eigenpy/decompositions/sparse/SimplicialLLT.hpp + include/eigenpy/decompositions/sparse/SimplicialLDLT.hpp + include/eigenpy/decompositions/sparse/SparseLU.hpp + include/eigenpy/decompositions/sparse/SparseQR.hpp include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp - include/eigenpy/decompositions/sparse/SparseSolverBase.hpp) + include/eigenpy/decompositions/sparse/SparseSolverBase.hpp + include/eigenpy/decompositions/sparse/LDLT.hpp + include/eigenpy/decompositions/sparse/LLT.hpp) if(BUILD_WITH_CHOLMOD_SUPPORT) list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS @@ -227,6 +236,16 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS ${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS} include/eigenpy/decompositions/decompositions.hpp include/eigenpy/decompositions/EigenSolver.hpp + include/eigenpy/decompositions/GeneralizedEigenSolver.hpp + include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp + include/eigenpy/decompositions/HessenbergDecomposition.hpp + include/eigenpy/decompositions/RealQZ.hpp + include/eigenpy/decompositions/Tridiagonalization.hpp + include/eigenpy/decompositions/RealSchur.hpp + include/eigenpy/decompositions/ComplexEigenSolver.hpp + include/eigenpy/decompositions/ComplexSchur.hpp + include/eigenpy/decompositions/FullPivLU.hpp + include/eigenpy/decompositions/PartialPivLU.hpp include/eigenpy/decompositions/PermutationMatrix.hpp include/eigenpy/decompositions/LDLT.hpp include/eigenpy/decompositions/LLT.hpp @@ -236,6 +255,9 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp include/eigenpy/decompositions/FullPivHouseholderQR.hpp include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp + include/eigenpy/decompositions/SVDBase.hpp + include/eigenpy/decompositions/BDCSVD.hpp + include/eigenpy/decompositions/JacobiSVD.hpp include/eigenpy/decompositions/minres.hpp) set(${PROJECT_NAME}_HEADERS @@ -300,17 +322,36 @@ list( # ---------------------------------------------------- # --- TARGETS ---------------------------------------- # ---------------------------------------------------- -set(${PROJECT_NAME}_SOLVERS_SOURCES src/solvers/preconditioners.cpp - src/solvers/solvers.cpp) +set(${PROJECT_NAME}_SOLVERS_SOURCES + src/solvers/preconditioners.cpp + src/solvers/solvers.cpp + src/solvers/minres.cpp + src/solvers/bicgstab.cpp + src/solvers/conjugate-gradient.cpp + src/solvers/least-squares-conjugate-gradient.cpp + src/solvers/incomplete-cholesky.cpp + src/solvers/incomplete-lut.cpp) set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES src/decompositions/decompositions.cpp src/decompositions/eigen-solver.cpp + src/decompositions/generalized-eigen-solver.cpp + src/decompositions/generalized-self-adjoint-eigen-solver.cpp + src/decompositions/complex-eigen-solver.cpp + src/decompositions/complex-schur.cpp src/decompositions/llt-solver.cpp src/decompositions/ldlt-solver.cpp - src/decompositions/minres-solver.cpp + src/decompositions/bdcsvd-solver.cpp + src/decompositions/jacobisvd-solver.cpp + src/decompositions/fullpivlu-solver.cpp + src/decompositions/hessenberg-decomposition.cpp + src/decompositions/real-qz.cpp + src/decompositions/tridiagonalization.cpp + src/decompositions/real-schur.cpp + src/decompositions/partialpivlu-solver.cpp + src/decompositions/sparse-lu-solver.cpp + src/decompositions/sparse-qr-solver.cpp src/decompositions/qr-solvers.cpp - src/decompositions/eigen-solver.cpp src/decompositions/self-adjoint-eigen-solver.cpp src/decompositions/permutation-matrix.cpp src/decompositions/simplicial-llt-solver.cpp diff --git a/README.md b/README.md index f5a5bae17..c1e9d8763 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,7 @@ EigenPy — Versatile and efficient Python bindings between Numpy and Eigen - full support Eigen::Ref avoiding memory allocation - full support of the Eigen::Tensor module - exposition of the Geometry module of Eigen for easy code prototyping -- standard matrix decomposion routines of Eigen such as the Cholesky decomposition (SVD and QR decompositions [can be added](#contributing)) +- standard matrix decomposion routines of Eigen such as the Cholesky, SVD and QR decompositions - full support of SWIG objects - full support of runtime declaration of Numpy scalar types - extended API to expose several STL types and some of their Boost equivalents: `optional` types, `std::pair`, maps, variants... @@ -100,6 +100,7 @@ The following people have been involved in the development of **EigenPy**: - [Loïc Estève](https://github.com/lesteve) (Inria): Conda integration - [Wilson Jallet](https://manifoldfr.github.io/) (Inria): core developer - [Joris Vaillant](https://github.com/jorisv) (Inria): core developer and manager of the project +- [Lucas Haubert](https://www.linkedin.com/in/lucas-haubert-b668a421a/) (Inria): core developer If you have taken part in the development of **EigenPy**, feel free to add your name and contribution here. diff --git a/flake.nix b/flake.nix index cd7e7d1f6..a0b48c61c 100644 --- a/flake.nix +++ b/flake.nix @@ -20,20 +20,29 @@ devShells.default = pkgs.mkShell { inputsFrom = [ self'.packages.default ]; }; packages = { default = self'.packages.eigenpy; - eigenpy = pkgs.python3Packages.eigenpy.overrideAttrs (_: { - src = pkgs.lib.fileset.toSource { - root = ./.; - fileset = pkgs.lib.fileset.unions [ - ./CMakeLists.txt - ./doc - ./include - ./package.xml - ./python - ./src - ./unittest - ]; - }; - }); + eigen = pkgs.eigen.overrideAttrs { + # Apply https://gitlab.com/libeigen/eigen/-/merge_requests/977 + postPatch = '' + substituteInPlace Eigen/src/SVD/BDCSVD.h \ + --replace-fail "if (l == 0) {" "if (i >= k && l == 0) {" + ''; + }; + eigenpy = + (pkgs.python3Packages.eigenpy.override { inherit (self'.packages) eigen; }).overrideAttrs + (_: { + src = pkgs.lib.fileset.toSource { + root = ./.; + fileset = pkgs.lib.fileset.unions [ + ./CMakeLists.txt + ./doc + ./include + ./package.xml + ./python + ./src + ./unittest + ]; + }; + }); }; }; }; diff --git a/include/eigenpy/decompositions/BDCSVD.hpp b/include/eigenpy/decompositions/BDCSVD.hpp new file mode 100644 index 000000000..9f9eb75ed --- /dev/null +++ b/include/eigenpy/decompositions/BDCSVD.hpp @@ -0,0 +1,99 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_bdcsvd_hpp__ +#define __eigenpy_decompositions_bdcsvd_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/eigen/EigenBase.hpp" +#include "eigenpy/decompositions/SVDBase.hpp" + +namespace eigenpy { + +template +struct BDCSVDVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef Eigen::BDCSVD Solver; + typedef typename MatrixType::Scalar Scalar; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init>( + bp::args("self", "rows", "cols", "computationOptions "), + "Default Constructor with memory preallocation. ")) + .def(bp::init>( + bp::args("self", "matrix", "computationOptions "), + "Constructor performing the decomposition of given matrix.")) + + .def("cols", &Solver::cols, bp::arg("self"), + "Returns the number of columns. ") + .def("compute", + (Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute, + bp::args("self", "matrix"), + "Method performing the decomposition of given matrix. Computes " + "Thin/Full " + "unitaries U/V if specified using the Options template parameter " + "or the class constructor. ", + bp::return_self<>()) + .def("compute", + (Solver & (Solver::*)(const MatrixType &matrix, + unsigned int computationOptions)) & + Solver::compute, + bp::args("self", "matrix", "computationOptions"), + "Method performing the decomposition of given matrix, as " + "specified by the computationOptions parameter. ", + bp::return_self<>()) + .def("rows", &Solver::rows, bp::arg("self"), + "Returns the number of rows. ") + .def("setSwitchSize", &Solver::setSwitchSize, bp::args("self", "s")) + + .def(SVDBaseVisitor()); + } + + static void expose() { + static const std::string classname = + "BDCSVD_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "Class Bidiagonal Divide and Conquer SVD.\n\n" + "This class first reduces the input matrix to bi-diagonal form using " + "class " + "UpperBidiagonalization, and then performs a divide-and-conquer " + "diagonalization. " + "Small blocks are diagonalized using class JacobiSVD. You can control " + "the " + "switching size with the setSwitchSize() method, default is 16. For " + "small matrice " + "(<16), it is thus preferable to directly use JacobiSVD. For larger " + "ones, BDCSVD " + "is highly recommended and can several order of magnitude faster.\n\n" + "Warming: this algorithm is unlikely to provide accurate result when " + "compiled with " + "unsafe math optimizations. For instance, this concerns Intel's " + "compiler (ICC), which " + "performs such optimization by default unless you compile with the " + "-fp-model precise " + "option. Likewise, the -ffast-math option of GCC or clang will " + "significantly degrade the " + "accuracy.", + bp::no_init) + .def(BDCSVDVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_bdcsvd_hpp__ diff --git a/include/eigenpy/decompositions/ComplexEigenSolver.hpp b/include/eigenpy/decompositions/ComplexEigenSolver.hpp new file mode 100644 index 000000000..691455418 --- /dev/null +++ b/include/eigenpy/decompositions/ComplexEigenSolver.hpp @@ -0,0 +1,86 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_complex_eigen_solver_hpp__ +#define __eigenpy_decompositions_complex_eigen_solver_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct ComplexEigenSolverVisitor : public boost::python::def_visitor< + ComplexEigenSolverVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::ComplexEigenSolver Solver; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>("Default constructor")) + .def(bp::init( + bp::arg("size"), "Default constructor with memory preallocation")) + .def(bp::init>( + bp::args("matrix", "compute_eigen_vectors"), + "Computes eigendecomposition of given matrix")) + + .def("eigenvalues", &Solver::eigenvalues, bp::arg("self"), + "Returns the eigenvalues of given matrix.", + bp::return_internal_reference<>()) + .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"), + "Returns the eigenvectors of given matrix.", + bp::return_internal_reference<>()) + + .def("compute", &ComplexEigenSolverVisitor::compute_proxy, + bp::args("self", "matrix"), + "Computes the eigendecomposition of given matrix.", + bp::return_self<>()) + .def("compute", + (Solver & + (Solver::*)(const Eigen::EigenBase& matrix, bool)) & + Solver::compute, + bp::args("self", "matrix", "compute_eigen_vectors"), + "Computes the eigendecomposition of given matrix.", + bp::return_self<>()) + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("getMaxIterations", &Solver::getMaxIterations, bp::arg("self"), + "Returns the maximum number of iterations.") + .def("setMaxIterations", &Solver::setMaxIterations, + bp::args("self", "max_iter"), + "Sets the maximum number of iterations allowed.", + bp::return_self<>()); + } + + static void expose() { + static const std::string classname = + "ComplexEigenSolver" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), bp::no_init) + .def(ComplexEigenSolverVisitor()) + .def(IdVisitor()); + } + + private: + template + static Solver& compute_proxy(Solver& self, + const Eigen::EigenBase& matrix) { + return self.compute(matrix); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_complex_eigen_solver_hpp__ diff --git a/include/eigenpy/decompositions/ComplexSchur.hpp b/include/eigenpy/decompositions/ComplexSchur.hpp new file mode 100644 index 000000000..bac2ea1bd --- /dev/null +++ b/include/eigenpy/decompositions/ComplexSchur.hpp @@ -0,0 +1,90 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_complex_schur_hpp__ +#define __eigenpy_decompositions_complex_schur_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct ComplexSchurVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::ComplexSchur Solver; + + template + void visit(PyClass& cl) const { + cl.def(bp::init(bp::arg("size"), "Default constructor")) + .def(bp::init>( + bp::args("matrix", "computeU"), "Computes Schur of given matrix")) + + .def("matrixU", &Solver::matrixU, bp::arg("self"), + "Returns the unitary matrix in the Schur decomposition. ", + bp::return_value_policy()) + .def("matrixT", &Solver::matrixT, bp::arg("self"), + "Returns the triangular matrix in the Schur decomposition. ", + bp::return_value_policy()) + + .def("compute", &ComplexSchurVisitor::compute_proxy, + bp::args("self", "matrix"), "Computes the Schur of given matrix.", + bp::return_self<>()) + .def("compute", + (Solver & + (Solver::*)(const Eigen::EigenBase& matrix, bool)) & + Solver::compute, + bp::args("self", "matrix", "computeU"), + "Computes the Schur of given matrix.", bp::return_self<>()) + + .def("computeFromHessenberg", + (Solver & (Solver::*)(const Eigen::EigenBase& matrixH, + const Eigen::EigenBase& matrixQ, + bool)) & + Solver::computeFromHessenberg, + bp::args("self", "matrixH", "matrixQ", "computeU"), + "Compute Schur decomposition from a given Hessenberg matrix. ", + bp::return_self<>()) + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("getMaxIterations", &Solver::getMaxIterations, bp::arg("self"), + "Returns the maximum number of iterations.") + .def("setMaxIterations", &Solver::setMaxIterations, + bp::args("self", "max_iter"), + "Sets the maximum number of iterations allowed.", + bp::return_self<>()); + } + + static void expose() { + static const std::string classname = + "ComplexSchur" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), bp::no_init) + .def(ComplexSchurVisitor()) + .def(IdVisitor()); + } + + private: + template + static Solver& compute_proxy(Solver& self, + const Eigen::EigenBase& matrix) { + return self.compute(matrix); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_complex_schur_hpp__ diff --git a/include/eigenpy/decompositions/FullPivLU.hpp b/include/eigenpy/decompositions/FullPivLU.hpp new file mode 100644 index 000000000..1a6c3788c --- /dev/null +++ b/include/eigenpy/decompositions/FullPivLU.hpp @@ -0,0 +1,192 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_fullpivlu_hpp__ +#define __eigenpy_decompositions_fullpivlu_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/eigen/EigenBase.hpp" + +namespace eigenpy { + +template +struct FullPivLUSolverVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix + VectorXs; + typedef Eigen::Matrix + MatrixXs; + typedef Eigen::FullPivLU Solver; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init( + bp::args("self", "rows", "cols"), + "Default constructor with memory preallocation")) + .def(bp::init( + bp::args("self", "matrix"), + "Constructs a LU factorization from a given matrix.")) + + .def(EigenBaseVisitor()) + + .def( + "compute", + (Solver & (Solver::*)(const Eigen::EigenBase &matrix)) & + Solver::compute, + bp::args("self", "matrix"), + "Computes the LU decomposition of the given matrix.", + bp::return_self<>()) + + .def("determinant", &Solver::determinant, bp::arg("self"), + "Returns the determinant of the matrix of which *this is the LU " + "decomposition.") + .def("dimensionOfKernel", &Solver::dimensionOfKernel, bp::arg("self"), + "Returns the dimension of the kernel of the matrix of which *this " + "is the LU decomposition.") + .def( + "image", + +[](Solver &self, const MatrixType &mat) -> MatrixType { + return self.image(mat); + }, + bp::args("self", "originalMatrix"), + "Returns the image of the matrix, also called its column-space. " + "The columns of the returned matrix will form a basis of the " + "image (column-space).") + .def( + "inverse", + +[](Solver &self) -> MatrixType { return self.inverse(); }, + bp::arg("self"), + "Returns the inverse of the matrix of which *this is the LU " + "decomposition.") + + .def("isInjective", &Solver::isInjective, bp::arg("self")) + .def("isInvertible", &Solver::isInvertible, bp::arg("self")) + .def("isSurjective", &Solver::isSurjective, bp::arg("self")) + + .def( + "kernel", +[](Solver &self) -> MatrixType { return self.kernel(); }, + bp::arg("self"), + "Returns the kernel of the matrix, also called its null-space. " + "The columns of the returned matrix will form a basis of the " + "kernel.") + + .def("matrixLU", &Solver::matrixLU, bp::arg("self"), + "Returns the LU decomposition matrix.", + bp::return_internal_reference<>()) + + .def("maxPivot", &Solver::maxPivot, bp::arg("self")) + .def("nonzeroPivots", &Solver::nonzeroPivots, bp::arg("self")) + + .def("permutationP", &Solver::permutationP, bp::arg("self"), + "Returns the permutation P.", + bp::return_value_policy()) + .def("permutationQ", &Solver::permutationQ, bp::arg("self"), + "Returns the permutation Q.", + bp::return_value_policy()) + + .def("rank", &Solver::rank, bp::arg("self")) + + .def("rcond", &Solver::rcond, bp::arg("self"), + "Returns an estimate of the reciprocal condition number of the " + "matrix.") + .def("reconstructedMatrix", &Solver::reconstructedMatrix, + bp::arg("self"), + "Returns the matrix represented by the decomposition, i.e., it " + "returns the product: P-1LUQ-1. This function is provided for " + "debug " + "purpose.") + + .def("setThreshold", + (Solver & (Solver::*)(const RealScalar &)) & Solver::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 LU 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( + "setThreshold", + +[](Solver &self) -> Solver & { + return self.setThreshold(Eigen::Default); + }, + bp::arg("self"), + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + bp::return_self<>()) + + .def("solve", &solve, bp::args("self", "b"), + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def("solve", &solve, 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.") + + .def("threshold", &Solver::threshold, bp::arg("self"), + "Returns the threshold that will be used by certain methods such " + "as rank()."); + } + + static void expose() { + static const std::string classname = + "FullPivLU" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "LU decomposition of a matrix with complete pivoting, and related " + "features.\n\n" + "This class represents a LU decomposition of any matrix, with complete " + "pivoting: the matrix A is decomposed as A=P−1LUQ−1 where L is " + "unit-lower-triangular, U is upper-triangular, and P and Q are " + "permutation matrices. This is a rank-revealing LU decomposition. " + "The eigenvalues (diagonal coefficients) of U are sorted in such a " + "way that any zeros are at the end.\n\n" + "This decomposition provides the generic approach to solving systems " + "of " + "linear equations, computing the rank, invertibility, inverse, kernel, " + "and determinant. \n\n" + "This LU decomposition is very stable and well tested with large " + "matrices. " + "However there are use cases where the SVD decomposition is inherently " + "more " + "stable and/or flexible. For example, when computing the kernel of a " + "matrix, " + "working with the SVD allows to select the smallest singular values of " + "the matrix, something that the LU decomposition doesn't see.", + bp::no_init) + .def(IdVisitor()) + .def(FullPivLUSolverVisitor()); + } + + private: + template + static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) { + return self.solve(vec); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_fullpivlu_hpp__ diff --git a/include/eigenpy/decompositions/GeneralizedEigenSolver.hpp b/include/eigenpy/decompositions/GeneralizedEigenSolver.hpp new file mode 100644 index 000000000..255470985 --- /dev/null +++ b/include/eigenpy/decompositions/GeneralizedEigenSolver.hpp @@ -0,0 +1,92 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_generalized_eigen_solver_hpp__ +#define __eigenpy_decompositions_generalized_eigen_solver_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct GeneralizedEigenSolverVisitor + : public boost::python::def_visitor< + GeneralizedEigenSolverVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::GeneralizedEigenSolver Solver; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>("Default constructor")) + .def(bp::init( + bp::arg("size"), "Default constructor with memory preallocation. ")) + .def(bp::init>( + bp::args("A", "B", "computeEigenVectors"), + "Computes the generalized eigendecomposition of given matrix " + "pair. ")) + + .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"), + "Returns an expression of the computed generalized eigenvectors. ") + .def( + "eigenvalues", + +[](const Solver& c) { return c.eigenvalues().eval(); }, + "Returns the computed generalized eigenvalues.") + + .def("alphas", &Solver::alphas, bp::arg("self"), + "Returns the vectors containing the alpha values. ") + .def("betas", &Solver::betas, bp::arg("self"), + "Returns the vectors containing the beta values. ") + + .def("compute", + &GeneralizedEigenSolverVisitor::compute_proxy, + bp::args("self", "A", "B"), + "Computes generalized eigendecomposition of given matrix. ", + bp::return_self<>()) + .def("compute", + (Solver & + (Solver::*)(const MatrixType& A, const MatrixType& B, bool)) & + Solver::compute, + bp::args("self", "A", "B", "computeEigenvectors"), + "Computes generalized eigendecomposition of given matrix. .", + bp::return_self<>()) + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("setMaxIterations", &Solver::setMaxIterations, + bp::args("self", "max_iter"), + "Sets the maximum number of iterations allowed.", + bp::return_self<>()); + } + + static void expose() { + static const std::string classname = + "GeneralizedEigenSolver" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), bp::no_init) + .def(GeneralizedEigenSolverVisitor()) + .def(IdVisitor()); + } + + private: + template + static Solver& compute_proxy(Solver& self, const MatrixType& A, + const MatrixType& B) { + return self.compute(A, B); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_generalized_eigen_solver_hpp__ diff --git a/include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp b/include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp new file mode 100644 index 000000000..7d23c7ac6 --- /dev/null +++ b/include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp @@ -0,0 +1,75 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__ +#define __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp" + +namespace eigenpy { + +template +struct GeneralizedSelfAdjointEigenSolverVisitor + : public boost::python::def_visitor< + GeneralizedSelfAdjointEigenSolverVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::GeneralizedSelfAdjointEigenSolver Solver; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>("Default constructor")) + .def(bp::init( + bp::arg("size"), "Default constructor with memory preallocation. ")) + .def(bp::init>( + bp::args("matA", "matB", "options"), + "Constructor: Computes generalized eigendecomposition of given " + "matrix pencil. ")) + + .def("compute", + &GeneralizedSelfAdjointEigenSolverVisitor::compute_proxy< + MatrixType>, + bp::args("self", "A", "B"), + "Computes generalized eigendecomposition of given matrix pencil. ", + bp::return_self<>()) + .def("compute", + (Solver & + (Solver::*)(const MatrixType& A, const MatrixType& B, int)) & + Solver::compute, + bp::args("self", "A", "B", "options"), + "Computes generalized eigendecomposition of given matrix pencil.", + bp::return_self<>()); + } + + static void expose() { + static const std::string classname = + "GeneralizedSelfAdjointEigenSolver" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_>>( + name.c_str(), bp::no_init) + .def(GeneralizedSelfAdjointEigenSolverVisitor()) + .def(IdVisitor()); + } + + private: + template + static Solver& compute_proxy(Solver& self, const MatrixType& A, + const MatrixType& B) { + return self.compute(A, B); + } +}; + +} // namespace eigenpy + +#endif // ifndef + // __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__ diff --git a/include/eigenpy/decompositions/HessenbergDecomposition.hpp b/include/eigenpy/decompositions/HessenbergDecomposition.hpp new file mode 100644 index 000000000..a547a3fb7 --- /dev/null +++ b/include/eigenpy/decompositions/HessenbergDecomposition.hpp @@ -0,0 +1,76 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_hessenberg_decomposition_hpp__ +#define __eigenpy_decompositions_hessenberg_decomposition_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct HessenbergDecompositionVisitor + : public boost::python::def_visitor< + HessenbergDecompositionVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::HessenbergDecomposition Solver; + + template + void visit(PyClass& cl) const { + cl + .def(bp::init( + bp::arg("size"), + "Default constructor; the decomposition will be computed later. ")) + .def(bp::init( + bp::arg("matrix"), + "Constructor; computes Hessenberg decomposition of given matrix. ")) + + .def( + "compute", + (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & + Solver::compute, + bp::args("self", "A"), + "Computes Hessenberg decomposition of given matrix. ", + bp::return_self<>()) + + .def("householderCoefficients", &Solver::householderCoefficients, + bp::arg("self"), "Returns the Householder coefficients. ", + bp::return_value_policy()) + + .def( + "matrixQ", + +[](const Solver& c) -> MatrixType { return c.matrixQ(); }, + "Reconstructs the orthogonal matrix Q in the decomposition.") + .def( + "matrixH", + +[](const Solver& c) -> MatrixType { return c.matrixH(); }, + "Constructs the Hessenberg matrix H in the decomposition.") + + .def("packedMatrix", &Solver::packedMatrix, bp::arg("self"), + "Returns the internal representation of the decomposition. ", + bp::return_value_policy()); + } + + static void expose() { + static const std::string classname = + "HessenbergDecomposition" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), bp::no_init) + .def(HessenbergDecompositionVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_hessenberg_decomposition_hpp__ diff --git a/include/eigenpy/decompositions/JacobiSVD.hpp b/include/eigenpy/decompositions/JacobiSVD.hpp new file mode 100644 index 000000000..1ea45ed1d --- /dev/null +++ b/include/eigenpy/decompositions/JacobiSVD.hpp @@ -0,0 +1,105 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_jacobisvd_hpp__ +#define __eigenpy_decompositions_jacobisvd_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/eigen/EigenBase.hpp" +#include "eigenpy/decompositions/SVDBase.hpp" + +namespace eigenpy { + +template +struct JacobiSVDVisitor + : public boost::python::def_visitor> { + typedef typename JacobiSVD::MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init>( + bp::args("self", "rows", "cols", "computationOptions "), + "Default Constructor with memory preallocation.")) + .def(bp::init>( + bp::args("self", "matrix", "computationOptions "), + "Constructor performing the decomposition of given matrix.")) + + .def("cols", &JacobiSVD::cols, bp::arg("self"), + "Returns the number of columns. ") + .def("compute", + (JacobiSVD & (JacobiSVD::*)(const MatrixType &matrix)) & + JacobiSVD::compute, + bp::args("self", "matrix"), + "Method performing the decomposition of given matrix. Computes " + "Thin/Full " + "unitaries U/V if specified using the Options template parameter " + "or the class constructor. ", + bp::return_self<>()) + .def("compute", + (JacobiSVD & (JacobiSVD::*)(const MatrixType &matrix, + unsigned int computationOptions)) & + JacobiSVD::compute, + bp::args("self", "matrix", "computationOptions"), + "Method performing the decomposition of given matrix, as " + "specified by the computationOptions parameter. ", + bp::return_self<>()) + .def("rows", &JacobiSVD::rows, bp::arg("self"), + "Returns the number of rows.") + + .def(SVDBaseVisitor()); + } + + static void expose() { + static const std::string classname = + "JacobiSVD_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "Two-sided Jacobi SVD decomposition of a rectangular matrix. \n\n" + "SVD decomposition consists in decomposing any n-by-p matrix A as a " + "product " + "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" + "eal positive matrix which is zero outside of its main diagonal; the " + "diagonal " + "entries of S are known as the singular values of A and the columns of " + "U and V " + "are known as the left and right singular vectors of A respectively. " + "\n\n" + "Singular values are always sorted in decreasing order. \n\n " + "This JacobiSVD decomposition computes only the singular values by " + "default. " + "If you want U or V, you need to ask for them explicitly. \n\n" + "You can ask for only thin U or V to be computed, meaning the " + "following. " + "In case of a rectangular n-by-p matrix, letting m be the smaller " + "value among " + "n and p, there are only m singular vectors; the remaining columns of " + "U and V " + "do not correspond to actual singular vectors. Asking for thin U or V " + "means asking " + "for only their m first columns to be formed. So U is then a n-by-m " + "matrix, and V " + "is then a p-by-m matrix. Notice that thin U and V are all you need " + "for (least squares) " + "solving.", + bp::no_init) + .def(JacobiSVDVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_jacobisvd_hpp__ diff --git a/include/eigenpy/decompositions/PartialPivLU.hpp b/include/eigenpy/decompositions/PartialPivLU.hpp new file mode 100644 index 000000000..2c986a831 --- /dev/null +++ b/include/eigenpy/decompositions/PartialPivLU.hpp @@ -0,0 +1,138 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_partialpivlu_hpp__ +#define __eigenpy_decompositions_partialpivlu_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/eigen/EigenBase.hpp" + +namespace eigenpy { + +template +struct PartialPivLUSolverVisitor : public boost::python::def_visitor< + PartialPivLUSolverVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix + VectorXs; + typedef Eigen::Matrix + MatrixXs; + typedef Eigen::PartialPivLU Solver; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init( + bp::args("self", "size"), + "Default constructor with memory preallocation")) + .def(bp::init( + bp::args("self", "matrix"), + "Constructs a LU factorization from a given matrix.")) + + .def(EigenBaseVisitor()) + + .def("determinant", &Solver::determinant, bp::arg("self"), + "Returns the determinant of the matrix of which *this is the LU " + "decomposition.") + .def( + "compute", + (Solver & (Solver::*)(const Eigen::EigenBase &matrix)) & + Solver::compute, + bp::args("self", "matrix"), + "Computes the LU factorization of given matrix.", + bp::return_self<>()) + .def( + "inverse", + +[](const Solver &self) -> MatrixType { return self.inverse(); }, + bp::arg("self"), + "Returns the inverse of the matrix of which *this is the LU " + "decomposition.") + .def("matrixLU", &Solver::matrixLU, bp::arg("self"), + "Returns the LU decomposition matrix.", + bp::return_internal_reference<>()) + + .def("permutationP", &Solver::permutationP, bp::arg("self"), + "Returns the permutation P.", + bp::return_value_policy()) + + .def("rcond", &Solver::rcond, bp::arg("self"), + "Returns an estimate of the reciprocal condition number of the " + "matrix.") + .def("reconstructedMatrix", &Solver::reconstructedMatrix, + bp::arg("self"), + "Returns the matrix represented by the decomposition, i.e., it " + "returns the product: P-1LUQ-1. This function is provided for " + "debug " + "purpose.") + .def("solve", &solve, bp::args("self", "b"), + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def("solve", &solve, 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 = + "FullPivLU" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "LU decomposition of a matrix with partial pivoting, " + "and related features. \n\n" + "This class represents a LU decomposition of a square " + "invertible matrix, " + "with partial pivoting: the matrix A is decomposed as A " + "= PLU where L is " + "unit-lower-triangular, U is upper-triangular, and P is " + "a permutation matrix.\n\n" + "Typically, partial pivoting LU decomposition is only " + "considered numerically " + "stable for square invertible matrices. Thus LAPACK's " + "dgesv and dgesvx require " + "the matrix to be square and invertible. The present " + "class does the same. It " + "will assert that the matrix is square, but it won't " + "(actually it can't) check " + "that the matrix is invertible: it is your task to " + "check that you only use this " + "decomposition on invertible matrices. \n\n" + "The guaranteed safe alternative, working for all matrices, " + "is the full pivoting LU decomposition, provided by class " + "FullPivLU. \n\n" + "This is not a rank-revealing LU decomposition. Many features " + "are intentionally absent from this class, such as " + "rank computation. If you need these features, use class " + "FullPivLU. \n\n" + "This LU decomposition is suitable to invert invertible " + "matrices. It is what MatrixBase::inverse() uses in the " + "general case. On the other hand, it is not suitable to " + "determine whether a given matrix is invertible. \n\n" + "The data of the LU decomposition can be directly accessed " + "through the methods matrixLU(), permutationP().", + bp::no_init) + .def(IdVisitor()) + .def(PartialPivLUSolverVisitor()); + } + + private: + template + static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) { + return self.solve(vec); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_partialpivlu_hpp__ diff --git a/include/eigenpy/decompositions/RealQZ.hpp b/include/eigenpy/decompositions/RealQZ.hpp new file mode 100644 index 000000000..76b45dbc2 --- /dev/null +++ b/include/eigenpy/decompositions/RealQZ.hpp @@ -0,0 +1,90 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_generalized_real_qz_hpp__ +#define __eigenpy_decompositions_generalized_real_qz_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct RealQZVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::RealQZ Solver; + + template + void visit(PyClass& cl) const { + cl.def( + bp::init(bp::arg("size"), "Default constructor. ")) + .def(bp::init>( + bp::args("A", "B", "computeQZ"), + "Constructor; computes real QZ decomposition of given matrices. ")) + + .def("compute", &RealQZVisitor::compute_proxy, + bp::args("self", "A", "B"), + "Computes QZ decomposition of given matrix. ", + bp::return_self<>()) + .def("compute", + (Solver & + (Solver::*)(const MatrixType& A, const MatrixType& B, bool)) & + Solver::compute, + bp::args("self", "A", "B", "computeEigenvectors"), + "Computes QZ decomposition of given matrix. ", bp::return_self<>()) + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("matrixQ", &Solver::matrixQ, bp::arg("self"), + "Returns matrix Q in the QZ decomposition. ", + bp::return_value_policy()) + .def("matrixS", &Solver::matrixS, bp::arg("self"), + "Returns matrix S in the QZ decomposition. ", + bp::return_value_policy()) + .def("matrixT", &Solver::matrixT, bp::arg("self"), + "Returns matrix T in the QZ decomposition. ", + bp::return_value_policy()) + .def("matrixZ", &Solver::matrixZ, bp::arg("self"), + "Returns matrix Z in the QZ decomposition. ", + bp::return_value_policy()) + + .def("iterations", &Solver::iterations, bp::arg("self"), + "Returns number of performed QR-like iterations. ") + .def("setMaxIterations", &Solver::setMaxIterations, + bp::args("self", "max_iter"), + "Sets the maximum number of iterations allowed.", + bp::return_self<>()); + } + + static void expose() { + static const std::string classname = + "RealQZVisitor" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), bp::no_init) + .def(RealQZVisitor()) + .def(IdVisitor()); + } + + private: + template + static Solver& compute_proxy(Solver& self, const MatrixType& A, + const MatrixType& B) { + return self.compute(A, B); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_generalized_real_qz_hpp__ diff --git a/include/eigenpy/decompositions/RealSchur.hpp b/include/eigenpy/decompositions/RealSchur.hpp new file mode 100644 index 000000000..ba0de9e59 --- /dev/null +++ b/include/eigenpy/decompositions/RealSchur.hpp @@ -0,0 +1,92 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_generalized_real_schur_hpp__ +#define __eigenpy_decompositions_generalized_real_schur_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct RealSchurVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::RealSchur Solver; + + template + void visit(PyClass& cl) const { + cl.def( + bp::init(bp::arg("size"), "Default constructor. ")) + .def(bp::init>( + bp::args("matrix", "computeU"), + "Constructor; computes real Schur decomposition of given matrix. ")) + + .def("compute", &RealSchurVisitor::compute_proxy, + bp::args("self", "matrix"), + "Computes Schur decomposition of given matrix. ", + bp::return_self<>()) + .def("compute", + (Solver & + (Solver::*)(const Eigen::EigenBase& matrix, bool)) & + Solver::compute, + bp::args("self", "matrix", "computeEigenvectors"), + "Computes Schur decomposition of given matrix. ", + bp::return_self<>()) + + .def("computeFromHessenberg", + (Solver & (Solver::*)(const MatrixType& matrixH, + const MatrixType& matrixQ, bool)) & + Solver::computeFromHessenberg, + bp::args("self", "matrixH", "matrixQ", "computeU"), + "Compute Schur decomposition from a given Hessenberg matrix. ", + bp::return_self<>()) + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("matrixT", &Solver::matrixT, bp::arg("self"), + "Returns the quasi-triangular matrix in the Schur decomposition.", + bp::return_value_policy()) + .def("matrixU", &Solver::matrixU, bp::arg("self"), + "Returns the orthogonal matrix in the Schur decomposition. ", + bp::return_value_policy()) + + .def("setMaxIterations", &Solver::setMaxIterations, + bp::args("self", "max_iter"), + "Sets the maximum number of iterations allowed.", + bp::return_self<>()) + .def("getMaxIterations", &Solver::getMaxIterations, bp::arg("self"), + "Returns the maximum number of iterations."); + } + + static void expose() { + static const std::string classname = + "RealSchurVisitor" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), bp::no_init) + .def(RealSchurVisitor()) + .def(IdVisitor()); + } + + private: + template + static Solver& compute_proxy(Solver& self, const MatrixType& A) { + return self.compute(A); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_generalized_real_schur_hpp__ diff --git a/include/eigenpy/decompositions/SVDBase.hpp b/include/eigenpy/decompositions/SVDBase.hpp new file mode 100644 index 000000000..5c9d0e391 --- /dev/null +++ b/include/eigenpy/decompositions/SVDBase.hpp @@ -0,0 +1,108 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_svdbase_hpp__ +#define __eigenpy_decompositions_svdbase_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/eigen/EigenBase.hpp" + +namespace eigenpy { + +template +struct SVDBaseVisitor + : public boost::python::def_visitor> { + typedef Derived Solver; + + typedef typename Derived::MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + + typedef Eigen::Matrix + VectorXs; + typedef Eigen::Matrix + MatrixXs; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + + .def("computeU", &Solver::computeU, bp::arg("self"), + "Returns true if U (full or thin) is asked for in this " + "SVD decomposition ") + .def("computeV", &Solver::computeV, bp::arg("self"), + "Returns true if V (full or thin) is asked for in this " + "SVD decomposition ") + + .def("info", &Solver::info, bp::arg("self"), + "Reports whether previous computation was successful. ") + + .def("matrixU", &matrixU, bp::arg("self"), "Returns the matrix U.") + .def("matrixV", &matrixV, bp::arg("self"), "Returns the matrix V.") + + .def("nonzeroSingularValues", &Solver::nonzeroSingularValues, + bp::arg("self"), + "Returns the number of singular values that are not exactly 0 ") + .def("rank", &Solver::rank, bp::arg("self"), + "the rank of the matrix of which *this is the SVD. ") + + .def("setThreshold", + (Solver & (Solver::*)(const RealScalar &)) & Solver::setThreshold, + bp::args("self", "threshold"), + "Allows to prescribe a threshold to be used by certain methods, " + "such as " + "rank() and solve(), which need to determine when singular values " + "are " + "to be considered nonzero. This is not used for the SVD " + "decomposition " + "itself.\n" + "\n" + "When it needs to get the threshold value, Eigen calls " + "threshold(). " + "The default is NumTraits::epsilon()", + bp::return_self<>()) + .def( + "setThreshold", + +[](Solver &self) -> Solver & { + return self.setThreshold(Eigen::Default); + }, + bp::arg("self"), + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + bp::return_self<>()) + + .def("singularValues", &Solver::singularValues, bp::arg("self"), + "Returns the vector of singular values.", + bp::return_value_policy()) + + .def("solve", &solve, bp::args("self", "b"), + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def("solve", &solve, 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.") + + .def("threshold", &Solver::threshold, bp::arg("self"), + "Returns the threshold that will be used by certain methods such " + "as rank(). "); + } + + private: + static MatrixXs matrixU(const Solver &self) { return self.matrixU(); } + static MatrixXs matrixV(const Solver &self) { return self.matrixV(); } + + template + static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) { + return self.solve(vec); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_svdbase_hpp__ diff --git a/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp b/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp index 221f3d434..987577d6f 100644 --- a/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp +++ b/include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp @@ -21,6 +21,7 @@ struct SelfAdjointEigenSolverVisitor typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef Eigen::SelfAdjointEigenSolver Solver; + typedef Eigen::Matrix VectorType; template void visit(PyClass& cl) const { @@ -32,12 +33,20 @@ struct SelfAdjointEigenSolverVisitor bp::args("self", "matrix", "options"), "Computes eigendecomposition of given matrix")) - .def("eigenvalues", &Solver::eigenvalues, bp::arg("self"), - "Returns the eigenvalues of given matrix.", - bp::return_internal_reference<>()) - .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"), - "Returns the eigenvectors of given matrix.", - bp::return_internal_reference<>()) + .def( + "eigenvalues", + +[](const Solver& c) -> const VectorType& { + return c.eigenvalues(); + }, + "Returns the eigenvalues of given matrix.", + bp::return_internal_reference<>()) + .def( + "eigenvectors", + +[](const Solver& c) -> const MatrixType& { + return c.eigenvectors(); + }, + "Returns the eigenvectors of given matrix.", + bp::return_internal_reference<>()) .def("compute", &SelfAdjointEigenSolverVisitor::compute_proxy, diff --git a/include/eigenpy/decompositions/Tridiagonalization.hpp b/include/eigenpy/decompositions/Tridiagonalization.hpp new file mode 100644 index 000000000..d24deb988 --- /dev/null +++ b/include/eigenpy/decompositions/Tridiagonalization.hpp @@ -0,0 +1,87 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_tridiagonalization_hpp__ +#define __eigenpy_decompositions_tridiagonalization_hpp__ + +#include +#include + +#include "eigenpy/eigen-to-python.hpp" +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct TridiagonalizationVisitor : public boost::python::def_visitor< + TridiagonalizationVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef Eigen::Tridiagonalization Solver; + typedef Eigen::VectorXd VectorType; + + template + void visit(PyClass &cl) const { + cl.def( + bp::init(bp::arg("size"), "Default constructor. ")) + .def(bp::init(bp::arg("matrix"), + "Constructor; computes tridiagonal " + "decomposition of given matrix. ")) + + .def( + "compute", + (Solver & (Solver::*)(const Eigen::EigenBase &matrix)) & + Solver::compute, + bp::args("self", "matrix"), + "Computes tridiagonal decomposition of given matrix. ", + bp::return_self<>()) + + .def("householderCoefficients", &Solver::householderCoefficients, + bp::arg("self"), "Returns the Householder coefficients. ") + .def("packedMatrix", &Solver::packedMatrix, bp::arg("self"), + "Returns the internal representation of the decomposition. ", + bp::return_value_policy()) + + .def( + "matrixQ", + +[](const Solver &c) -> MatrixType { return c.matrixQ(); }, + "Returns the unitary matrix Q in the decomposition.") + .def( + "matrixT", + +[](const Solver &c) -> MatrixType { return c.matrixT(); }, + "Returns an expression of the tridiagonal matrix T in the " + "decomposition.") + + .def( + "diagonal", + +[](const Solver &c) -> VectorType { return c.diagonal(); }, + bp::arg("self"), + "Returns the diagonal of the tridiagonal matrix T in the " + "decomposition. ") + + .def( + "subDiagonal", + +[](const Solver &c) -> VectorType { return c.subDiagonal(); }, + bp::arg("self"), + "Returns the subdiagonal of the tridiagonal matrix T in the " + "decomposition."); + } + + static void expose() { + static const std::string classname = + "TridiagonalizationVisitor" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_(name.c_str(), bp::no_init) + .def(TridiagonalizationVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_tridiagonalization_hpp__ diff --git a/include/eigenpy/decompositions/minres.hpp b/include/eigenpy/decompositions/minres.hpp index fa53870e9..c5bb0dae0 100644 --- a/include/eigenpy/decompositions/minres.hpp +++ b/include/eigenpy/decompositions/minres.hpp @@ -1,181 +1,16 @@ /* - * Copyright 2021 INRIA + * Copyright 2025 INRIA */ #ifndef __eigenpy_decompositions_minres_hpp__ #define __eigenpy_decompositions_minres_hpp__ -#include -#include -#include +#include "eigenpy/fwd.hpp" -#include "eigenpy/eigenpy.hpp" -#include "eigenpy/utils/scalar-name.hpp" +// clang-format off +EIGENPY_PRAGMA_DEPRECATED_HEADER(eigenpy/decompositions/minres.hpp, eigenpy/solvers/MINRES.hpp) +// clang-format on -namespace eigenpy { -template -struct IterativeSolverBaseVisitor - : public boost::python::def_visitor> { - typedef _Solver Solver; - typedef typename Solver::MatrixType MatrixType; - typedef typename Solver::Preconditioner Preconditioner; - typedef typename Solver::Scalar Scalar; - typedef typename Solver::RealScalar RealScalar; - - typedef Eigen::Matrix - MatrixXs; - - template - void visit(PyClass& cl) const { - cl.def("analyzePattern", - (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & - Solver::analyzePattern, - bp::args("self", "A"), - "Initializes the iterative solver for the sparsity pattern of the " - "matrix A for further solving Ax=b problems.", - bp::return_self<>()) - .def( - "factorize", - (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & - Solver::factorize, - bp::args("self", "A"), - "Initializes the iterative solver with the numerical values of the " - "matrix A for further solving Ax=b problems.", - bp::return_self<>()) - .def( - "compute", - (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & - Solver::compute, - bp::args("self", "A"), - "Initializes the iterative solver with the matrix A for further " - "solving Ax=b problems.", - bp::return_self<>()) - - .def("rows", &Solver::rows, bp::arg("self"), - "Returns the number of rows.") - .def("cols", &Solver::cols, bp::arg("self"), - "Returns the number of columns.") - .def("tolerance", &Solver::tolerance, bp::arg("self"), - "Returns the tolerance threshold used by the stopping criteria.") - .def("setTolerance", &Solver::setTolerance, - bp::args("self", "tolerance"), - "Sets the tolerance threshold used by the stopping criteria.\n" - "This value is used as an upper bound to the relative residual " - "error: |Ax-b|/|b|.\n" - "The default value is the machine precision given by " - "NumTraits::epsilon().", - bp::return_self<>()) - .def("preconditioner", - (Preconditioner & (Solver::*)()) & Solver::preconditioner, - bp::arg("self"), - "Returns a read-write reference to the preconditioner for custom " - "configuration.", - bp::return_internal_reference<>()) - - .def("maxIterations", &Solver::maxIterations, bp::arg("self"), - "Returns the max number of iterations.\n" - "It is either the value setted by setMaxIterations or, by " - "default, twice the number of columns of the matrix.") - .def("setMaxIterations", &Solver::setMaxIterations, - bp::args("self", "max_iterations"), - "Sets the max number of iterations.\n" - "Default is twice the number of columns of the matrix.", - bp::return_self<>()) - - .def( - "iterations", &Solver::iterations, bp::arg("self"), - "Returns the number of iterations performed during the last solve.") - .def("error", &Solver::error, bp::arg("self"), - "Returns the tolerance error reached during the last solve.\n" - "It is a close approximation of the true relative residual error " - "|Ax-b|/|b|.") - .def("info", &Solver::error, bp::arg("info"), - "Returns Success if the iterations converged, and NoConvergence " - "otherwise.") - - .def("solveWithGuess", &solveWithGuess, - bp::args("self", "b", "x0"), - "Returns the solution x of A x = b using the current " - "decomposition of A and x0 as an initial solution.") - - .def( - "solve", &solve, 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 or vector."); - } - - private: - template - static MatrixOrVector1 solveWithGuess(const Solver& self, - const MatrixOrVector1& b, - const MatrixOrVector2& guess) { - return self.solveWithGuess(b, guess); - } - - template - static MatrixOrVector solve(const Solver& self, - const MatrixOrVector& mat_or_vec) { - MatrixOrVector res = self.solve(mat_or_vec); - return res; - } -}; - -template -struct MINRESSolverVisitor - : public boost::python::def_visitor> { - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef Eigen::Matrix - VectorXs; - typedef Eigen::Matrix - MatrixXs; - typedef Eigen::MINRES Solver; - - template - void visit(PyClass& cl) const { - cl.def(bp::init<>(bp::arg("self"), "Default constructor")) - .def(bp::init( - bp::args("self", "matrix"), - "Initialize the solver with matrix A for further Ax=b solving.\n" - "This constructor is a shortcut for the default constructor " - "followed by a call to compute().")) - - .def(IterativeSolverBaseVisitor()); - } - - static void expose() { - static const std::string classname = - "MINRES" + scalar_name::shortname(); - expose(classname); - } - - static void expose(const std::string& name) { - bp::class_( - name.c_str(), - "A minimal residual solver for sparse symmetric problems.\n" - "This class allows to solve for A.x = b sparse linear problems using " - "the MINRES algorithm of Paige and Saunders (1975). The sparse matrix " - "A must be symmetric (possibly indefinite). The vectors x and b can be " - "either dense or sparse.\n" - "The maximal number of iterations and tolerance value can be " - "controlled via the setMaxIterations() and setTolerance() methods. The " - "defaults are the size of the problem for the maximal number of " - "iterations and NumTraits::epsilon() for the tolerance.\n", - bp::no_init) - .def(MINRESSolverVisitor()) - .def(IdVisitor()); - } - - private: - template - static MatrixOrVector solve(const Solver& self, const MatrixOrVector& vec) { - return self.solve(vec); - } -}; - -} // namespace eigenpy +#include "eigenpy/solvers/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 6c9c39773..07664e992 100644 --- a/include/eigenpy/decompositions/sparse/LDLT.hpp +++ b/include/eigenpy/decompositions/sparse/LDLT.hpp @@ -1,72 +1,16 @@ /* - * Copyright 2024 INRIA + * Copyright 2025 INRIA */ #ifndef __eigenpy_decompositions_sparse_ldlt_hpp__ #define __eigenpy_decompositions_sparse_ldlt_hpp__ -#include "eigenpy/eigenpy.hpp" -#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp" -#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/fwd.hpp" -namespace eigenpy { +// clang-format off +EIGENPY_PRAGMA_DEPRECATED_HEADER(eigenpy/decompositions/sparse/LDLT.hpp, eigenpy/decompositions/sparse/SimplicialLDLT.hpp) +// clang-format on -template > -struct SimplicialLDLTVisitor - : public boost::python::def_visitor< - SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering>> { - typedef SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor; - typedef _MatrixType MatrixType; - - typedef Eigen::SimplicialLDLT Solver; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef Eigen::Matrix - DenseVectorXs; - typedef Eigen::Matrix - DenseMatrixXs; - - template - void visit(PyClass &cl) const { - cl.def(bp::init<>(bp::arg("self"), "Default constructor")) - .def(bp::init(bp::args("self", "matrix"), - "Constructs and performs the LDLT " - "factorization from a given matrix.")) - - .def("vectorD", &vectorD, bp::arg("self"), - "Returns the diagonal vector D.") - .def(SimplicialCholeskyVisitor()); - } - - static void expose() { - static const std::string classname = - "SimplicialLDLT_" + scalar_name::shortname(); - expose(classname); - } - - static void expose(const std::string &name) { - bp::class_( - name.c_str(), - "A direct sparse LDLT Cholesky factorizations.\n\n" - "This class provides a LDL^T Cholesky factorizations of sparse " - "matrices that are selfadjoint and positive definite." - "The factorization allows for solving A.X = B where X and B can be " - "either dense or sparse.\n\n" - "In order to reduce the fill-in, a symmetric permutation P is applied " - "prior to the factorization such that the factorized matrix is P A " - "P^-1.", - bp::no_init) - .def(SimplicialLDLTVisitor()) - .def(IdVisitor()); - } - - private: - static DenseVectorXs vectorD(const Solver &self) { return self.vectorD(); } -}; - -} // namespace eigenpy +#include "eigenpy/decompositions/sparse/SimplicialLDLT.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 259f796c1..795d18e5c 100644 --- a/include/eigenpy/decompositions/sparse/LLT.hpp +++ b/include/eigenpy/decompositions/sparse/LLT.hpp @@ -1,67 +1,16 @@ /* - * Copyright 2024 INRIA + * Copyright 2025 INRIA */ #ifndef __eigenpy_decompositions_sparse_llt_hpp__ #define __eigenpy_decompositions_sparse_llt_hpp__ -#include "eigenpy/eigenpy.hpp" -#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp" -#include "eigenpy/utils/scalar-name.hpp" +#include "eigenpy/fwd.hpp" -namespace eigenpy { +// clang-format off +EIGENPY_PRAGMA_DEPRECATED_HEADER(eigenpy/decompositions/sparse/LLT.hpp, eigenpy/decompositions/sparse/SimplicialLLT.hpp) +// clang-format on -template > -struct SimplicialLLTVisitor - : public boost::python::def_visitor< - SimplicialLLTVisitor<_MatrixType, _UpLo, _Ordering>> { - typedef SimplicialLLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor; - typedef _MatrixType MatrixType; - - typedef Eigen::SimplicialLLT Solver; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef Eigen::Matrix - DenseVectorXs; - typedef Eigen::Matrix - DenseMatrixXs; - - template - void visit(PyClass &cl) const { - cl.def(bp::init<>(bp::arg("self"), "Default constructor")) - .def(bp::init(bp::args("self", "matrix"), - "Constructs and performs the LLT " - "factorization from a given matrix.")) - - .def(SimplicialCholeskyVisitor()); - } - - static void expose() { - static const std::string classname = - "SimplicialLLT_" + scalar_name::shortname(); - expose(classname); - } - - static void expose(const std::string &name) { - bp::class_( - name.c_str(), - "A direct sparse LLT Cholesky factorizations.\n\n" - "This class provides a LL^T Cholesky factorizations of sparse matrices " - "that are selfadjoint and positive definite." - "The factorization allows for solving A.X = B where X and B can be " - "either dense or sparse.\n\n" - "In order to reduce the fill-in, a symmetric permutation P is applied " - "prior to the factorization such that the factorized matrix is P A " - "P^-1.", - bp::no_init) - .def(SimplicialLLTVisitor()) - .def(IdVisitor()); - } -}; - -} // namespace eigenpy +#include "eigenpy/decompositions/sparse/SimplicialLLT.hpp" #endif // ifndef __eigenpy_decompositions_sparse_llt_hpp__ diff --git a/include/eigenpy/decompositions/sparse/SimplicialLDLT.hpp b/include/eigenpy/decompositions/sparse/SimplicialLDLT.hpp new file mode 100644 index 000000000..3ba78b0fb --- /dev/null +++ b/include/eigenpy/decompositions/sparse/SimplicialLDLT.hpp @@ -0,0 +1,72 @@ +/* + * Copyright 2024 INRIA + */ + +#ifndef __eigenpy_decompositions_sparse_simplicial_ldlt_hpp__ +#define __eigenpy_decompositions_sparse_simplicial_ldlt_hpp__ + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template > +struct SimplicialLDLTVisitor + : public boost::python::def_visitor< + SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering>> { + typedef SimplicialLDLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor; + typedef _MatrixType MatrixType; + + typedef Eigen::SimplicialLDLT Solver; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix + DenseVectorXs; + typedef Eigen::Matrix + DenseMatrixXs; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init(bp::args("self", "matrix"), + "Constructs and performs the LDLT " + "factorization from a given matrix.")) + + .def("vectorD", &vectorD, bp::arg("self"), + "Returns the diagonal vector D.") + .def(SimplicialCholeskyVisitor()); + } + + static void expose() { + static const std::string classname = + "SimplicialLDLT_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "A direct sparse LDLT Cholesky factorizations.\n\n" + "This class provides a LDL^T Cholesky factorizations of sparse " + "matrices that are selfadjoint and positive definite." + "The factorization allows for solving A.X = B where X and B can be " + "either dense or sparse.\n\n" + "In order to reduce the fill-in, a symmetric permutation P is applied " + "prior to the factorization such that the factorized matrix is P A " + "P^-1.", + bp::no_init) + .def(SimplicialLDLTVisitor()) + .def(IdVisitor()); + } + + private: + static DenseVectorXs vectorD(const Solver &self) { return self.vectorD(); } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_sparse_simplicial_ldlt_hpp__ diff --git a/include/eigenpy/decompositions/sparse/SimplicialLLT.hpp b/include/eigenpy/decompositions/sparse/SimplicialLLT.hpp new file mode 100644 index 000000000..50581e9cd --- /dev/null +++ b/include/eigenpy/decompositions/sparse/SimplicialLLT.hpp @@ -0,0 +1,67 @@ +/* + * Copyright 2024 INRIA + */ + +#ifndef __eigenpy_decompositions_sparse_simplicial_llt_hpp__ +#define __eigenpy_decompositions_sparse_simplicial_llt_hpp__ + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template > +struct SimplicialLLTVisitor + : public boost::python::def_visitor< + SimplicialLLTVisitor<_MatrixType, _UpLo, _Ordering>> { + typedef SimplicialLLTVisitor<_MatrixType, _UpLo, _Ordering> Visitor; + typedef _MatrixType MatrixType; + + typedef Eigen::SimplicialLLT Solver; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix + DenseVectorXs; + typedef Eigen::Matrix + DenseMatrixXs; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init(bp::args("self", "matrix"), + "Constructs and performs the LLT " + "factorization from a given matrix.")) + + .def(SimplicialCholeskyVisitor()); + } + + static void expose() { + static const std::string classname = + "SimplicialLLT_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "A direct sparse LLT Cholesky factorizations.\n\n" + "This class provides a LL^T Cholesky factorizations of sparse matrices " + "that are selfadjoint and positive definite." + "The factorization allows for solving A.X = B where X and B can be " + "either dense or sparse.\n\n" + "In order to reduce the fill-in, a symmetric permutation P is applied " + "prior to the factorization such that the factorized matrix is P A " + "P^-1.", + bp::no_init) + .def(SimplicialLLTVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_sparse_simplicial_llt_hpp__ diff --git a/include/eigenpy/decompositions/sparse/SparseLU.hpp b/include/eigenpy/decompositions/sparse/SparseLU.hpp new file mode 100644 index 000000000..045b3f44e --- /dev/null +++ b/include/eigenpy/decompositions/sparse/SparseLU.hpp @@ -0,0 +1,217 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_sparse_lu_hpp__ +#define __eigenpy_decompositions_sparse_lu_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/decompositions/sparse/SparseSolverBase.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct SparseLUMatrixLReturnTypeVisitor + : public boost::python::def_visitor< + SparseLUMatrixLReturnTypeVisitor> { + typedef Eigen::SparseLUMatrixLReturnType LType; + typedef typename MappedSupernodalType::Scalar Scalar; + typedef Eigen::Matrix VectorXs; + typedef Eigen::Matrix + MatrixXs; + + template + static void solveInPlace(const LType &self, + Eigen::Ref mat_vec) { + self.solveInPlace(mat_vec); + } + + template + void visit(PyClass &cl) const { + cl.def(bp::init(bp::args("self", "mapL"))) + + .def("rows", <ype::rows) + .def("cols", <ype::cols) + + .def("solveInPlace", &solveInPlace, bp::args("self", "X")) + .def("solveInPlace", &solveInPlace, bp::args("self", "x")); + } + + static void expose(const std::string &name) { + bp::class_(name.c_str(), "Eigen SparseLUMatrixLReturnType", + bp::no_init) + .def(SparseLUMatrixLReturnTypeVisitor()) + .def(IdVisitor()); + } +}; + +template +struct SparseLUMatrixUReturnTypeVisitor + : public boost::python::def_visitor< + SparseLUMatrixUReturnTypeVisitor> { + typedef Eigen::SparseLUMatrixUReturnType UType; + typedef typename MatrixLType::Scalar Scalar; + typedef Eigen::Matrix VectorXs; + typedef Eigen::Matrix + MatrixXs; + + template + static void solveInPlace(const UType &self, + Eigen::Ref mat_vec) { + self.solveInPlace(mat_vec); + } + + template + void visit(PyClass &cl) const { + cl.def(bp::init(bp::args("self", "mapL", "mapU"))) + + .def("rows", &UType::rows) + .def("cols", &UType::cols) + + .def("solveInPlace", &solveInPlace, bp::args("self", "X")) + .def("solveInPlace", &solveInPlace, bp::args("self", "x")); + } + + static void expose(const std::string &name) { + bp::class_(name.c_str(), "Eigen SparseLUMatrixUReturnType", + bp::no_init) + .def(SparseLUMatrixUReturnTypeVisitor()) + .def(IdVisitor()); + } +}; + +template > +struct SparseLUVisitor : public boost::python::def_visitor< + SparseLUVisitor<_MatrixType, _Ordering>> { + typedef SparseLUVisitor<_MatrixType, _Ordering> Visitor; + typedef _MatrixType MatrixType; + + typedef Eigen::SparseLU Solver; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + + typedef typename Solver::SCMatrix SCMatrix; + typedef typename MatrixType::StorageIndex StorageIndex; + typedef Eigen::MappedSparseMatrix + MappedSparseMatrix; + typedef Eigen::SparseLUMatrixLReturnType LType; + typedef Eigen::SparseLUMatrixUReturnType UType; + + template + void visit(PyClass &cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init(bp::args("self", "matrix"), + "Constructs and performs the LU " + "factorization from a given matrix.")) + + .def("determinant", &Solver::determinant, bp::arg("self"), + "Returns the determinant of the matrix.") + .def("signDeterminant", &Solver::signDeterminant, bp::arg("self"), + "A number representing the sign of the determinant. ") + .def("absDeterminant", &Solver::absDeterminant, bp::arg("self"), + "Returns the absolute value of the determinant of the matrix of " + "which *this is the QR decomposition.") + .def("logAbsDeterminant", &Solver::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") + + .def("rows", &Solver::rows, bp::arg("self"), + "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, bp::arg("self"), + "Returns the number of cols of the matrix.") + + .def("nnzL", &Solver::nnzL, bp::arg("self"), + "The number of non zero elements in L") + .def("nnzU", &Solver::nnzU, bp::arg("self"), + "The number of non zero elements in U") + + .def( + "matrixL", + +[](const Solver &self) -> LType { return self.matrixL(); }, + "Returns an expression of the matrix L.") + .def( + "matrixU", + +[](const Solver &self) -> UType { return self.matrixU(); }, + "Returns an expression of the matrix U.") + + .def("colsPermutation", &Solver::colsPermutation, bp::arg("self"), + "Returns a reference to the column matrix permutation PTc such " + "that Pr A PTc = LU.", + bp::return_value_policy()) + .def("rowsPermutation", &Solver::rowsPermutation, bp::arg("self"), + "Returns a reference to the row matrix permutation Pr such that " + "Pr A PTc = LU", + bp::return_value_policy()) + + .def("compute", &Solver::compute, bp::args("self", "matrix"), + "Compute the symbolic and numeric factorization of the input " + "sparse matrix. " + "The input matrix should be in column-major storage. ") + + .def("analyzePattern", &Solver::analyzePattern, + bp::args("self", "matrix"), + "Compute the column permutation to minimize the fill-in.") + .def("factorize", &Solver::factorize, bp::args("self", "matrix"), + "Performs a numeric decomposition of a given matrix.\n" + "The given matrix must has the same sparcity than the matrix on " + "which the symbolic decomposition has been performed.") + + .def("isSymmetric", &Solver::isSymmetric, bp::args("self", "sym"), + "Indicate that the pattern of the input matrix is symmetric. ") + + .def("setPivotThreshold", &Solver::setPivotThreshold, + bp::args("self", "thresh"), + "Set the threshold used for a diagonal entry to be an acceptable " + "pivot.") + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + .def("lastErrorMessage", &Solver::lastErrorMessage, bp::arg("self"), + "Returns a string describing the type of error. ") + + .def(SparseSolverBaseVisitor()); + } + + static void expose() { + static const std::string classname = + "SparseLU_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string &name) { + bp::class_( + name.c_str(), + "Sparse supernodal LU factorization for general matrices. \n\n" + "This class implements the supernodal LU factorization for general " + "matrices. " + "It uses the main techniques from the sequential SuperLU package " + "(http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles " + "transparently real " + "and complex arithmetic with single and double precision, depending on " + "the scalar " + "type of your input matrix. The code has been optimized to provide " + "BLAS-3 operations " + "during supernode-panel updates. It benefits directly from the " + "built-in high-performant " + "Eigen BLAS routines. Moreover, when the size of a supernode is very " + "small, the BLAS " + "calls are avoided to enable a better optimization from the compiler. " + "For best performance, " + "you should compile it with NDEBUG flag to avoid the numerous bounds " + "checking on vectors.", + bp::no_init) + .def(SparseLUVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_sparse_lu_hpp__ diff --git a/include/eigenpy/decompositions/sparse/SparseQR.hpp b/include/eigenpy/decompositions/sparse/SparseQR.hpp new file mode 100644 index 000000000..e78ca0025 --- /dev/null +++ b/include/eigenpy/decompositions/sparse/SparseQR.hpp @@ -0,0 +1,239 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_decompositions_sparse_qr_hpp__ +#define __eigenpy_decompositions_sparse_qr_hpp__ + +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/decompositions/sparse/SparseSolverBase.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct SparseQRMatrixQTransposeReturnTypeVisitor + : public boost::python::def_visitor< + SparseQRMatrixQTransposeReturnTypeVisitor> { + typedef typename SparseQRType::Scalar Scalar; + typedef Eigen::SparseQRMatrixQTransposeReturnType + QTransposeType; + typedef Eigen::VectorXd VectorXd; + typedef Eigen::MatrixXd MatrixXd; + + template + void visit(PyClass& cl) const { + cl.def(bp::init(bp::args("self", "qr"))) + .def( + "__matmul__", + +[](QTransposeType& self, const MatrixXd& other) -> MatrixXd { + return MatrixXd(self * other); + }, + bp::args("self", "other")) + + .def( + "__matmul__", + +[](QTransposeType& self, const VectorXd& other) -> VectorXd { + return VectorXd(self * other); + }, + bp::args("self", "other")); + } + + static void expose() { + static const std::string classname = "SparseQRMatrixQTransposeReturnType_" + + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_( + name.c_str(), "Eigen SparseQRMatrixQTransposeReturnType", bp::no_init) + .def(SparseQRMatrixQTransposeReturnTypeVisitor()) + .def(IdVisitor()); + } +}; + +template +struct SparseQRMatrixQReturnTypeVisitor + : public boost::python::def_visitor< + SparseQRMatrixQReturnTypeVisitor> { + typedef typename SparseQRType::Scalar Scalar; + typedef Eigen::SparseQRMatrixQTransposeReturnType + QTransposeType; + typedef Eigen::SparseQRMatrixQReturnType QType; + typedef typename SparseQRType::QRMatrixType QRMatrixType; + typedef Eigen::VectorXd VectorXd; + typedef Eigen::MatrixXd MatrixXd; + + template + void visit(PyClass& cl) const { + cl.def(bp::init(bp::args("self", "qr"))) + .def( + "__matmul__", + +[](QType& self, const MatrixXd& other) -> MatrixXd { + return MatrixXd(self * other); + }, + bp::args("self", "other")) + + .def( + "__matmul__", + +[](QType& self, const VectorXd& other) -> VectorXd { + return VectorXd(self * other); + }, + bp::args("self", "other")) + + .def("rows", &QType::rows) + .def("cols", &QType::cols) + + .def( + "adjoint", + +[](const QType& self) -> QTransposeType { return self.adjoint(); }) + + .def( + "transpose", + +[](const QType& self) -> QTransposeType { + return self.transpose(); + }) + + .def( + "toSparse", + +[](QType& self) -> QRMatrixType { + Eigen::Index m = self.rows(); + MatrixXd I = MatrixXd::Identity(m, m); + MatrixXd Q_dense = self * I; + return Q_dense.sparseView(); + }, + bp::arg("self"), + "Convert the Q matrix to a sparse matrix representation."); + } + + static void expose() { + static const std::string classname = + "SparseQRMatrixQReturnType_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), "Eigen SparseQRMatrixQReturnType", + bp::no_init) + .def(SparseQRMatrixQReturnTypeVisitor()) + .def(IdVisitor()); + } +}; + +template +struct SparseQRVisitor + : public boost::python::def_visitor> { + typedef typename SparseQRType::MatrixType MatrixType; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix + DenseVectorXs; + typedef Eigen::Matrix + DenseMatrixXs; + + typedef typename SparseQRType::QRMatrixType QRMatrixType; + typedef Eigen::SparseQRMatrixQReturnType QType; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init( + bp::args("self", "mat"), + "Construct a QR factorization of the matrix mat.")) + + .def("cols", &SparseQRType::cols, bp::arg("self"), + "Returns the number of columns of the represented matrix. ") + .def("rows", &SparseQRType::rows, bp::arg("self"), + "Returns the number of rows of the represented matrix. ") + + .def("compute", &SparseQRType::compute, bp::args("self", "matrix"), + "Compute the symbolic and numeric factorization of the input " + "sparse matrix. " + "The input matrix should be in column-major storage. ") + .def("analyzePattern", &SparseQRType::analyzePattern, + bp::args("self", "mat"), + "Compute the column permutation to minimize the fill-in.") + .def("factorize", &SparseQRType::factorize, bp::args("self", "matrix"), + "Performs a numeric decomposition of a given matrix.\n" + "The given matrix must has the same sparcity than the matrix on " + "which the symbolic decomposition has been performed.") + + .def("colsPermutation", &SparseQRType::colsPermutation, bp::arg("self"), + "Returns a reference to the column matrix permutation PTc such " + "that Pr A PTc = LU.", + bp::return_value_policy()) + + .def("info", &SparseQRType::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + .def("lastErrorMessage", &SparseQRType::lastErrorMessage, + bp::arg("self"), "Returns a string describing the type of error. ") + + .def("rank", &SparseQRType::rank, bp::arg("self"), + "Returns the number of non linearly dependent columns as " + "determined " + "by the pivoting threshold. ") + + .def( + "matrixQ", + +[](const SparseQRType& self) -> QType { return self.matrixQ(); }, + "Returns an expression of the matrix Q as products of sparse " + "Householder reflectors.") + .def( + "matrixR", + +[](const SparseQRType& self) -> const QRMatrixType& { + return self.matrixR(); + }, + "Returns a const reference to the \b sparse upper triangular " + "matrix " + "R of the QR factorization.", + bp::return_value_policy()) + + .def("setPivotThreshold", &SparseQRType::setPivotThreshold, + bp::args("self", "thresh"), + "Set the threshold used for a diagonal entry to be an acceptable " + "pivot.") + + .def(SparseSolverBaseVisitor()); + } + + static void expose() { + static const std::string classname = + "SparseQR_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_( + name.c_str(), + "Sparse left-looking QR factorization with numerical column pivoting. " + "This class implements a left-looking QR decomposition of sparse " + "matrices " + "with numerical column pivoting. When a column has a norm less than a " + "given " + "tolerance it is implicitly permuted to the end. The QR factorization " + "thus " + "obtained is given by A*P = Q*R where R is upper triangular or " + "trapezoidal. \n\n" + "P is the column permutation which is the product of the fill-reducing " + "and the " + "numerical permutations. \n\n" + "Q is the orthogonal matrix represented as products of Householder " + "reflectors. \n\n" + "R is the sparse triangular or trapezoidal matrix. The later occurs " + "when A is rank-deficient. \n\n", + bp::no_init) + .def(SparseQRVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decompositions_sparse_qr_hpp__ diff --git a/include/eigenpy/decompositions/sparse/accelerate/Accelerate.hpp b/include/eigenpy/decompositions/sparse/accelerate/Accelerate.hpp new file mode 100644 index 000000000..29cf9cdef --- /dev/null +++ b/include/eigenpy/decompositions/sparse/accelerate/Accelerate.hpp @@ -0,0 +1,73 @@ +/* + * Copyright 2024 INRIA + */ + +#ifndef __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__ +#define __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__ + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/eigen/EigenBase.hpp" +#include "eigenpy/decompositions/sparse/SparseSolverBase.hpp" + +#include + +namespace eigenpy { + +template +struct AccelerateImplVisitor : public boost::python::def_visitor< + AccelerateImplVisitor> { + typedef AccelerateDerived Solver; + + typedef typename AccelerateDerived::MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef MatrixType CholMatrixType; + typedef typename MatrixType::StorageIndex StorageIndex; + + template + void visit(PyClass &cl) const { + cl + + .def("analyzePattern", &Solver::analyzePattern, + bp::args("self", "matrix"), + "Performs a symbolic decomposition on the sparcity of matrix.\n" + "This function is particularly useful when solving for several " + "problems having the same structure.") + + .def(EigenBaseVisitor()) + .def(SparseSolverBaseVisitor()) + + .def("compute", + (Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute, + bp::args("self", "matrix"), + "Computes the sparse Cholesky decomposition of a given matrix.", + bp::return_self<>()) + + .def("factorize", &Solver::factorize, bp::args("self", "matrix"), + "Performs a numeric decomposition of a given matrix.\n" + "The given matrix must has the same sparcity than the matrix on " + "which the symbolic decomposition has been performed.\n" + "See also analyzePattern().") + + .def("info", &Solver::info, bp::arg("self"), + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("setOrder", &Solver::setOrder, bp::arg("self"), "Set order"); + } + + static void expose(const std::string &name, const std::string &doc = "") { + bp::class_(name.c_str(), doc.c_str(), + bp::no_init) + .def(AccelerateImplVisitor()) + + .def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init(bp::args("self", "matrix"), + "Constructs and performs the " + "factorization from a given matrix.")); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__ diff --git a/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp b/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp index 29cf9cdef..71ff88b05 100644 --- a/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp +++ b/include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp @@ -1,73 +1,18 @@ /* - * Copyright 2024 INRIA + * Copyright 2025 INRIA */ -#ifndef __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__ -#define __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__ +#ifndef __eigenpy_decompositions_sparse_accelerate_accelerate_hpp_deprecated__ +#define __eigenpy_decompositions_sparse_accelerate_accelerate_hpp_deprecated__ -#include "eigenpy/eigenpy.hpp" -#include "eigenpy/eigen/EigenBase.hpp" -#include "eigenpy/decompositions/sparse/SparseSolverBase.hpp" +#include "eigenpy/fwd.hpp" -#include +// clang-format off +EIGENPY_PRAGMA_DEPRECATED_HEADER(eigenpy/decompositions/sparse/accelerate/accelerate.hpp, eigenpy/decompositions/sparse/accelerate/Accelerate.hpp) +// clang-format on -namespace eigenpy { +#include "eigenpy/decompositions/sparse/accelerate/Accelerate.hpp" -template -struct AccelerateImplVisitor : public boost::python::def_visitor< - AccelerateImplVisitor> { - typedef AccelerateDerived Solver; +#endif - typedef typename AccelerateDerived::MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef MatrixType CholMatrixType; - typedef typename MatrixType::StorageIndex StorageIndex; - - template - void visit(PyClass &cl) const { - cl - - .def("analyzePattern", &Solver::analyzePattern, - bp::args("self", "matrix"), - "Performs a symbolic decomposition on the sparcity of matrix.\n" - "This function is particularly useful when solving for several " - "problems having the same structure.") - - .def(EigenBaseVisitor()) - .def(SparseSolverBaseVisitor()) - - .def("compute", - (Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute, - bp::args("self", "matrix"), - "Computes the sparse Cholesky decomposition of a given matrix.", - bp::return_self<>()) - - .def("factorize", &Solver::factorize, bp::args("self", "matrix"), - "Performs a numeric decomposition of a given matrix.\n" - "The given matrix must has the same sparcity than the matrix on " - "which the symbolic decomposition has been performed.\n" - "See also analyzePattern().") - - .def("info", &Solver::info, bp::arg("self"), - "NumericalIssue if the input contains INF or NaN values or " - "overflow occured. Returns Success otherwise.") - - .def("setOrder", &Solver::setOrder, bp::arg("self"), "Set order"); - } - - static void expose(const std::string &name, const std::string &doc = "") { - bp::class_(name.c_str(), doc.c_str(), - bp::no_init) - .def(AccelerateImplVisitor()) - - .def(bp::init<>(bp::arg("self"), "Default constructor")) - .def(bp::init(bp::args("self", "matrix"), - "Constructs and performs the " - "factorization from a given matrix.")); - } -}; - -} // namespace eigenpy - -#endif // ifndef __eigenpy_decomposition_sparse_accelerate_accelerate_hpp__ +// ifndef __eigenpy_decompositions_sparse_accelerate_accelerate_hpp_deprecated__ diff --git a/include/eigenpy/fwd.hpp b/include/eigenpy/fwd.hpp index c0fe4c75f..d95a81f2c 100644 --- a/include/eigenpy/fwd.hpp +++ b/include/eigenpy/fwd.hpp @@ -36,11 +36,11 @@ #if defined(EIGENPY_CLANG_COMPILER) || defined(EIGENPY_GCC_COMPILER) #define EIGENPY_PRAGMA(x) _Pragma(#x) #define EIGENPY_PRAGMA_MESSAGE(the_message) \ - EIGENPY_PRAGMA(GCC message the_message) + EIGENPY_PRAGMA(GCC message #the_message) #define EIGENPY_PRAGMA_WARNING(the_message) \ - EIGENPY_PRAGMA(GCC warning the_message) + EIGENPY_PRAGMA(GCC warning #the_message) #define EIGENPY_PRAGMA_DEPRECATED(the_message) \ - EIGENPY_PRAGMA_WARNING(Deprecated : the_message) + EIGENPY_PRAGMA_WARNING(Deprecated : #the_message) #define EIGENPY_PRAGMA_DEPRECATED_HEADER(old_header, new_header) \ EIGENPY_PRAGMA_WARNING( \ Deprecated header file : #old_header has been replaced \ diff --git a/include/eigenpy/solvers/BFGSPreconditioners.hpp b/include/eigenpy/solvers/BFGSPreconditioners.hpp index 40da2b584..2c591bb34 100644 --- a/include/eigenpy/solvers/BFGSPreconditioners.hpp +++ b/include/eigenpy/solvers/BFGSPreconditioners.hpp @@ -1,10 +1,10 @@ /* * Copyright 2017 CNRS - * Copyright 2024 Inria + * Copyright 2024-2025 INRIA */ -#ifndef __eigenpy_bfgs_preconditioners_hpp__ -#define __eigenpy_bfgs_preconditioners_hpp__ +#ifndef __eigenpy_solvers_bfgs_preconditioners_hpp__ +#define __eigenpy_solvers_bfgs_preconditioners_hpp__ #include @@ -65,4 +65,4 @@ struct LimitedBFGSPreconditionerBaseVisitor } // namespace eigenpy -#endif // ifndef __eigenpy_bfgs_preconditioners_hpp__ +#endif // ifndef __eigenpy_solvers_bfgs_preconditioners_hpp__ diff --git a/include/eigenpy/solvers/BasicPreconditioners.hpp b/include/eigenpy/solvers/BasicPreconditioners.hpp index 7e4cd7190..2f4a0a242 100644 --- a/include/eigenpy/solvers/BasicPreconditioners.hpp +++ b/include/eigenpy/solvers/BasicPreconditioners.hpp @@ -1,10 +1,10 @@ /* * Copyright 2017 CNRS - * Copyright 2024 Inria + * Copyright 2024-2025 INRIA */ -#ifndef __eigenpy_basic_preconditioners_hpp__ -#define __eigenpy_basic_preconditioners_hpp__ +#ifndef __eigenpy_solvers_basic_preconditioners_hpp__ +#define __eigenpy_solvers_basic_preconditioners_hpp__ #include @@ -115,4 +115,4 @@ struct IdentityPreconditionerVisitor } // namespace eigenpy -#endif // ifndef __eigenpy_basic_preconditioners_hpp__ +#endif // ifndef __eigenpy_solvers_basic_preconditioners_hpp__ diff --git a/include/eigenpy/solvers/BiCGSTAB.hpp b/include/eigenpy/solvers/BiCGSTAB.hpp new file mode 100644 index 000000000..41bd4bf3c --- /dev/null +++ b/include/eigenpy/solvers/BiCGSTAB.hpp @@ -0,0 +1,41 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_solvers_bicgstab_hpp__ +#define __eigenpy_solvers_bicgstab_hpp__ + +#include + +#include "eigenpy/fwd.hpp" +#include "eigenpy/solvers/IterativeSolverBase.hpp" + +namespace eigenpy { + +template +struct BiCGSTABVisitor + : public boost::python::def_visitor> { + typedef typename BiCGSTAB::MatrixType MatrixType; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>("Default constructor")) + .def(bp::init( + bp::arg("A"), + "Initialize the solver with matrix A for further || Ax - b || " + "solving.\n" + "This constructor is a shortcut for the default constructor " + "followed by a call to compute().")); + } + + static void expose(const std::string& name = "BiCGSTAB") { + bp::class_(name.c_str(), bp::no_init) + .def(IterativeSolverVisitor()) + .def(BiCGSTABVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_solvers_bicgstab_hpp__ diff --git a/include/eigenpy/solvers/ConjugateGradient.hpp b/include/eigenpy/solvers/ConjugateGradient.hpp index f2fdf8699..2c99b277a 100644 --- a/include/eigenpy/solvers/ConjugateGradient.hpp +++ b/include/eigenpy/solvers/ConjugateGradient.hpp @@ -1,9 +1,10 @@ /* * Copyright 2017 CNRS + * Copyright 2025 INRIA */ -#ifndef __eigenpy_conjugate_gradient_hpp__ -#define __eigenpy_conjugate_gradient_hpp__ +#ifndef __eigenpy_solvers_conjugate_gradient_hpp__ +#define __eigenpy_solvers_conjugate_gradient_hpp__ #include @@ -38,4 +39,4 @@ struct ConjugateGradientVisitor } // namespace eigenpy -#endif // ifndef __eigenpy_conjugate_gradient_hpp__ +#endif // ifndef __eigenpy_solvers_conjugate_gradient_hpp__ diff --git a/include/eigenpy/solvers/IncompleteCholesky.hpp b/include/eigenpy/solvers/IncompleteCholesky.hpp new file mode 100644 index 000000000..f0f1e30ca --- /dev/null +++ b/include/eigenpy/solvers/IncompleteCholesky.hpp @@ -0,0 +1,125 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_solvers_incomplete_cholesky_hpp__ +#define __eigenpy_solvers_incomplete_cholesky_hpp__ + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct IncompleteCholeskyVisitor : public boost::python::def_visitor< + IncompleteCholeskyVisitor<_MatrixType>> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::IncompleteCholesky Solver; + typedef Eigen::Matrix + DenseVectorXs; + typedef Eigen::Matrix + DenseMatrixXs; + + typedef Eigen::SparseMatrix FactorType; + typedef Eigen::Matrix VectorRx; + typedef Eigen::PermutationMatrix + PermutationType; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init(bp::args("self", "matrix"), + "Constructs and performs the LDLT " + "factorization from a given matrix.")) + + .def("rows", &Solver::rows, bp::arg("self"), + "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, bp::arg("self"), + "Returns the number of cols of the matrix.") + + .def("info", &Solver::info, bp::arg("self"), + "Reports whether previous computation was successful.") + + .def("setInitialShift", &Solver::setInitialShift, + bp::args("self", "shift"), "Set the initial shift parameter.") + + .def( + "analyzePattern", + +[](Solver& self, const MatrixType& amat) { + self.analyzePattern(amat); + }, + bp::arg("matrix")) + .def( + "factorize", + +[](Solver& self, const MatrixType& amat) { self.factorize(amat); }, + bp::arg("matrix")) + .def( + "compute", + +[](Solver& self, const MatrixType& amat) { self.compute(amat); }, + bp::arg("matrix")) + + .def( + "matrixL", + +[](const Solver& self) -> const FactorType& { + return self.matrixL(); + }, + bp::return_value_policy()) + .def( + "scalingS", + +[](const Solver& self) -> const VectorRx& { + return self.scalingS(); + }, + bp::return_value_policy()) + .def( + "permutationP", + +[](const Solver& self) -> const PermutationType& { + return self.permutationP(); + }, + bp::return_value_policy()) + + .def( + "solve", + +[](const Solver& self, const Eigen::Ref& b) + -> DenseVectorXs { return self.solve(b); }, + bp::arg("b"), + "Returns the solution x of A x = b using the current decomposition " + "of A, where b is a right hand side vector.") + .def( + "solve", + +[](const Solver& self, const Eigen::Ref& B) + -> DenseMatrixXs { return self.solve(B); }, + bp::arg("b"), + "Returns the solution X of A X = B using the current decomposition " + "of A where B is a right hand side matrix.") + .def( + "solve", + +[](const Solver& self, const MatrixType& B) -> MatrixType { + DenseMatrixXs B_dense = DenseMatrixXs(B); + DenseMatrixXs X_dense = self.solve(B_dense); + return MatrixType(X_dense.sparseView()); + }, + bp::arg("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 = + "IncompleteCholesky_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_(name.c_str(), "Incomplete Cholesky.", + bp::no_init) + .def(IncompleteCholeskyVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_solvers_incomplete_cholesky_hpp__ diff --git a/include/eigenpy/solvers/IncompleteLUT.hpp b/include/eigenpy/solvers/IncompleteLUT.hpp new file mode 100644 index 000000000..04c58ff74 --- /dev/null +++ b/include/eigenpy/solvers/IncompleteLUT.hpp @@ -0,0 +1,108 @@ +/* + * Copyright 2025 INRIA + */ + +#ifndef __eigenpy_solvers_incomplete_lut_hpp__ +#define __eigenpy_solvers_incomplete_lut_hpp__ + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { + +template +struct IncompleteLUTVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::IncompleteLUT Solver; + typedef Eigen::Matrix + DenseVectorXs; + typedef Eigen::Matrix + DenseMatrixXs; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init(bp::args("self", "matrix"), + "Constructs and performs the LDLT " + "factorization from a given matrix.")) + .def(bp::init( + (bp::arg("matrix"), + bp::arg("droptol") = Eigen::NumTraits::dummy_precision(), + bp::arg("fillfactor") = 10), + "Constructs an incomplete LU factorization from a given matrix.")) + + .def("rows", &Solver::rows, bp::arg("self"), + "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, bp::arg("self"), + "Returns the number of cols of the matrix.") + + .def("info", &Solver::info, bp::arg("self"), + "Reports whether previous computation was successful.") + + .def( + "analyzePattern", + +[](Solver& self, const MatrixType& amat) { + self.analyzePattern(amat); + }, + bp::arg("matrix")) + .def( + "factorize", + +[](Solver& self, const MatrixType& amat) { self.factorize(amat); }, + bp::arg("matrix")) + .def( + "compute", + +[](Solver& self, const MatrixType& amat) { self.compute(amat); }, + bp::arg("matrix")) + + .def("setDroptol", &Solver::setDroptol, bp::arg("self")) + .def("setFillfactor", &Solver::setFillfactor, bp::arg("self")) + + .def( + "solve", + +[](const Solver& self, const Eigen::Ref& b) + -> DenseVectorXs { return self.solve(b); }, + bp::arg("b"), + "Returns the solution x of A x = b using the current decomposition " + "of A, where b is a right hand side vector.") + .def( + "solve", + +[](const Solver& self, const Eigen::Ref& B) + -> DenseMatrixXs { return self.solve(B); }, + bp::arg("b"), + "Returns the solution X of A X = B using the current decomposition " + "of A where B is a right hand side matrix.") + .def( + "solve", + +[](const Solver& self, const MatrixType& B) -> MatrixType { + DenseMatrixXs B_dense = DenseMatrixXs(B); + DenseMatrixXs X_dense = self.solve(B_dense); + return MatrixType(X_dense.sparseView()); + }, + bp::arg("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 = + "IncompleteLUT_" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_( + name.c_str(), + "Incomplete LU factorization with dual-threshold strategy.", + bp::no_init) + .def(IncompleteLUTVisitor()) + .def(IdVisitor()); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_solvers_incomplete_lut_hpp__ diff --git a/include/eigenpy/solvers/IterativeSolverBase.hpp b/include/eigenpy/solvers/IterativeSolverBase.hpp index 670bde3d2..ebab5763a 100644 --- a/include/eigenpy/solvers/IterativeSolverBase.hpp +++ b/include/eigenpy/solvers/IterativeSolverBase.hpp @@ -1,9 +1,10 @@ /* * Copyright 2017 CNRS + * Copyright 2025 INRIA */ -#ifndef __eigenpy_iterative_solver_base_hpp__ -#define __eigenpy_iterative_solver_base_hpp__ +#ifndef __eigenpy_solvers_iterative_solver_base_hpp__ +#define __eigenpy_solvers_iterative_solver_base_hpp__ #include "eigenpy/fwd.hpp" #include "eigenpy/solvers/SparseSolverBase.hpp" @@ -103,4 +104,4 @@ struct IterativeSolverVisitor : public boost::python::def_visitor< } // namespace eigenpy -#endif // ifndef __eigenpy_iterative_solver_base_hpp__ +#endif // ifndef __eigenpy_solvers_iterative_solver_base_hpp__ diff --git a/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp b/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp index 043d3d7fb..2f3aae644 100644 --- a/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp +++ b/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp @@ -1,9 +1,10 @@ /* * Copyright 2017-2018 CNRS + * Copyright 2025 INRIA */ -#ifndef __eigenpy_least_square_conjugate_gradient_hpp__ -#define __eigenpy_least_square_conjugate_gradient_hpp__ +#ifndef __eigenpy_solvers_least_square_conjugate_gradient_hpp__ +#define __eigenpy_solvers_least_square_conjugate_gradient_hpp__ #include @@ -16,7 +17,7 @@ template struct LeastSquaresConjugateGradientVisitor : public boost::python::def_visitor< LeastSquaresConjugateGradientVisitor> { - typedef Eigen::MatrixXd MatrixType; + typedef typename LeastSquaresConjugateGradient::MatrixType MatrixType; template void visit(PyClass& cl) const { @@ -29,9 +30,10 @@ struct LeastSquaresConjugateGradientVisitor "followed by a call to compute().")); } - static void expose() { - bp::class_( - "LeastSquaresConjugateGradient", bp::no_init) + static void expose( + const std::string& name = "LeastSquaresConjugateGradient") { + bp::class_(name.c_str(), + bp::no_init) .def(IterativeSolverVisitor()) .def(LeastSquaresConjugateGradientVisitor< LeastSquaresConjugateGradient>()) @@ -41,4 +43,4 @@ struct LeastSquaresConjugateGradientVisitor } // namespace eigenpy -#endif // ifndef __eigenpy_least_square_conjugate_gradient_hpp__ +#endif // ifndef __eigenpy_solvers_least_square_conjugate_gradient_hpp__ diff --git a/include/eigenpy/solvers/MINRES.hpp b/include/eigenpy/solvers/MINRES.hpp new file mode 100644 index 000000000..c21fcb66d --- /dev/null +++ b/include/eigenpy/solvers/MINRES.hpp @@ -0,0 +1,181 @@ +/* + * Copyright 2021-2025 INRIA + */ + +#ifndef __eigenpy_solvers_minres_hpp__ +#define __eigenpy_solvers_minres_hpp__ + +#include +#include +#include + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy { +template +struct IterativeSolverBaseVisitor + : public boost::python::def_visitor> { + typedef _Solver Solver; + typedef typename Solver::MatrixType MatrixType; + typedef typename Solver::Preconditioner Preconditioner; + typedef typename Solver::Scalar Scalar; + typedef typename Solver::RealScalar RealScalar; + + typedef Eigen::Matrix + MatrixXs; + + template + void visit(PyClass& cl) const { + cl.def("analyzePattern", + (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & + Solver::analyzePattern, + bp::args("self", "A"), + "Initializes the iterative solver for the sparsity pattern of the " + "matrix A for further solving Ax=b problems.", + bp::return_self<>()) + .def( + "factorize", + (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & + Solver::factorize, + bp::args("self", "A"), + "Initializes the iterative solver with the numerical values of the " + "matrix A for further solving Ax=b problems.", + bp::return_self<>()) + .def( + "compute", + (Solver & (Solver::*)(const Eigen::EigenBase& matrix)) & + Solver::compute, + bp::args("self", "A"), + "Initializes the iterative solver with the matrix A for further " + "solving Ax=b problems.", + bp::return_self<>()) + + .def("rows", &Solver::rows, bp::arg("self"), + "Returns the number of rows.") + .def("cols", &Solver::cols, bp::arg("self"), + "Returns the number of columns.") + .def("tolerance", &Solver::tolerance, bp::arg("self"), + "Returns the tolerance threshold used by the stopping criteria.") + .def("setTolerance", &Solver::setTolerance, + bp::args("self", "tolerance"), + "Sets the tolerance threshold used by the stopping criteria.\n" + "This value is used as an upper bound to the relative residual " + "error: |Ax-b|/|b|.\n" + "The default value is the machine precision given by " + "NumTraits::epsilon().", + bp::return_self<>()) + .def("preconditioner", + (Preconditioner & (Solver::*)()) & Solver::preconditioner, + bp::arg("self"), + "Returns a read-write reference to the preconditioner for custom " + "configuration.", + bp::return_internal_reference<>()) + + .def("maxIterations", &Solver::maxIterations, bp::arg("self"), + "Returns the max number of iterations.\n" + "It is either the value setted by setMaxIterations or, by " + "default, twice the number of columns of the matrix.") + .def("setMaxIterations", &Solver::setMaxIterations, + bp::args("self", "max_iterations"), + "Sets the max number of iterations.\n" + "Default is twice the number of columns of the matrix.", + bp::return_self<>()) + + .def( + "iterations", &Solver::iterations, bp::arg("self"), + "Returns the number of iterations performed during the last solve.") + .def("error", &Solver::error, bp::arg("self"), + "Returns the tolerance error reached during the last solve.\n" + "It is a close approximation of the true relative residual error " + "|Ax-b|/|b|.") + .def("info", &Solver::error, bp::arg("info"), + "Returns Success if the iterations converged, and NoConvergence " + "otherwise.") + + .def("solveWithGuess", &solveWithGuess, + bp::args("self", "b", "x0"), + "Returns the solution x of A x = b using the current " + "decomposition of A and x0 as an initial solution.") + + .def( + "solve", &solve, 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 or vector."); + } + + private: + template + static MatrixOrVector1 solveWithGuess(const Solver& self, + const MatrixOrVector1& b, + const MatrixOrVector2& guess) { + return self.solveWithGuess(b, guess); + } + + template + static MatrixOrVector solve(const Solver& self, + const MatrixOrVector& mat_or_vec) { + MatrixOrVector res = self.solve(mat_or_vec); + return res; + } +}; + +template +struct MINRESSolverVisitor + : public boost::python::def_visitor> { + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix + VectorXs; + typedef Eigen::Matrix + MatrixXs; + typedef Eigen::MINRES Solver; + + template + void visit(PyClass& cl) const { + cl.def(bp::init<>(bp::arg("self"), "Default constructor")) + .def(bp::init( + bp::args("self", "matrix"), + "Initialize the solver with matrix A for further Ax=b solving.\n" + "This constructor is a shortcut for the default constructor " + "followed by a call to compute().")) + + .def(IterativeSolverBaseVisitor()); + } + + static void expose() { + static const std::string classname = + "MINRES" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string& name) { + bp::class_( + name.c_str(), + "A minimal residual solver for sparse symmetric problems.\n" + "This class allows to solve for A.x = b sparse linear problems using " + "the MINRES algorithm of Paige and Saunders (1975). The sparse matrix " + "A must be symmetric (possibly indefinite). The vectors x and b can be " + "either dense or sparse.\n" + "The maximal number of iterations and tolerance value can be " + "controlled via the setMaxIterations() and setTolerance() methods. The " + "defaults are the size of the problem for the maximal number of " + "iterations and NumTraits::epsilon() for the tolerance.\n", + bp::no_init) + .def(MINRESSolverVisitor()) + .def(IdVisitor()); + } + + private: + template + static MatrixOrVector solve(const Solver& self, const MatrixOrVector& vec) { + return self.solve(vec); + } +}; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_solvers_minres_hpp__ diff --git a/include/eigenpy/solvers/SparseSolverBase.hpp b/include/eigenpy/solvers/SparseSolverBase.hpp index f4f2bdad1..795d98eef 100644 --- a/include/eigenpy/solvers/SparseSolverBase.hpp +++ b/include/eigenpy/solvers/SparseSolverBase.hpp @@ -1,9 +1,10 @@ /* * Copyright 2017 CNRS + * Copyright 2025 INRIA */ -#ifndef __eigenpy_sparse_solver_base_hpp__ -#define __eigenpy_sparse_solver_base_hpp__ +#ifndef __eigenpy_solvers_sparse_solver_base_hpp__ +#define __eigenpy_solvers_sparse_solver_base_hpp__ #include "eigenpy/fwd.hpp" @@ -29,4 +30,4 @@ struct SparseSolverVisitor } // namespace eigenpy -#endif // ifndef __eigenpy_sparse_solver_base_hpp__ +#endif // ifndef __eigenpy_solvers_sparse_solver_base_hpp__ diff --git a/include/eigenpy/solvers/preconditioners.hpp b/include/eigenpy/solvers/preconditioners.hpp index 3e26638ba..c537e2779 100644 --- a/include/eigenpy/solvers/preconditioners.hpp +++ b/include/eigenpy/solvers/preconditioners.hpp @@ -1,9 +1,10 @@ /* * Copyright 2017 CNRS + * Copyright 2025 INRIA */ -#ifndef __eigenpy_preconditioners_hpp__ -#define __eigenpy_preconditioners_hpp__ +#ifndef __eigenpy_solvers_preconditioners_hpp__ +#define __eigenpy_solvers_preconditioners_hpp__ #include "eigenpy/config.hpp" @@ -13,4 +14,4 @@ void EIGENPY_DLLAPI exposePreconditioners(); } // namespace eigenpy -#endif // define __eigenpy_preconditioners_hpp__ +#endif // define __eigenpy_solvers_preconditioners_hpp__ diff --git a/include/eigenpy/solvers/solvers.hpp b/include/eigenpy/solvers/solvers.hpp index af1a375c3..16e1fddce 100644 --- a/include/eigenpy/solvers/solvers.hpp +++ b/include/eigenpy/solvers/solvers.hpp @@ -1,9 +1,9 @@ /* - * Copyright 2017-2020 CNRS INRIA + * Copyright 2017-2025 CNRS INRIA */ -#ifndef __eigenpy_solvers_hpp__ -#define __eigenpy_solvers_hpp__ +#ifndef __eigenpy_solvers_solvers_hpp__ +#define __eigenpy_solvers_solvers_hpp__ #include "eigenpy/config.hpp" @@ -14,4 +14,4 @@ void EIGENPY_DLLAPI exposeSolvers(); } // namespace eigenpy -#endif // define __eigenpy_solvers_hpp__ +#endif // define __eigenpy_solvers_solvers_hpp__ diff --git a/src/decompositions/accelerate.cpp b/src/decompositions/accelerate.cpp index 8f2f05c16..ae688a2e0 100644 --- a/src/decompositions/accelerate.cpp +++ b/src/decompositions/accelerate.cpp @@ -5,7 +5,7 @@ #include "eigenpy/fwd.hpp" #include "eigenpy/decompositions/decompositions.hpp" -#include "eigenpy/decompositions/sparse/accelerate/accelerate.hpp" +#include "eigenpy/decompositions/sparse/accelerate/Accelerate.hpp" namespace eigenpy { diff --git a/src/decompositions/bdcsvd-solver.cpp b/src/decompositions/bdcsvd-solver.cpp new file mode 100644 index 000000000..2fda1b3fb --- /dev/null +++ b/src/decompositions/bdcsvd-solver.cpp @@ -0,0 +1,12 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/BDCSVD.hpp" + +namespace eigenpy { +void exposeBDCSVDSolver() { + using namespace Eigen; + BDCSVDVisitor::expose("BDCSVD"); +} +} // namespace eigenpy diff --git a/src/decompositions/complex-eigen-solver.cpp b/src/decompositions/complex-eigen-solver.cpp new file mode 100644 index 000000000..99df394ad --- /dev/null +++ b/src/decompositions/complex-eigen-solver.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/ComplexEigenSolver.hpp" + +namespace eigenpy { +void exposeComplexEigenSolver() { + using namespace Eigen; + ComplexEigenSolverVisitor::expose("ComplexEigenSolver"); +} +} // namespace eigenpy diff --git a/src/decompositions/complex-schur.cpp b/src/decompositions/complex-schur.cpp new file mode 100644 index 000000000..555e345c2 --- /dev/null +++ b/src/decompositions/complex-schur.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/ComplexSchur.hpp" + +namespace eigenpy { +void exposeComplexSchur() { + using namespace Eigen; + ComplexSchurVisitor::expose("ComplexSchur"); +} +} // namespace eigenpy diff --git a/src/decompositions/decompositions.cpp b/src/decompositions/decompositions.cpp index 77ae22f29..85fb81061 100644 --- a/src/decompositions/decompositions.cpp +++ b/src/decompositions/decompositions.cpp @@ -3,30 +3,60 @@ */ #include "eigenpy/decompositions/decompositions.hpp" +#include "eigenpy/solvers/MINRES.hpp" #include "eigenpy/fwd.hpp" namespace eigenpy { void exposeEigenSolver(); +void exposeGeneralizedEigenSolver(); void exposeSelfAdjointEigenSolver(); +void exposeGeneralizedSelfAdjointEigenSolver(); +void exposeHessenbergDecomposition(); +void exposeRealQZ(); +void exposeRealSchur(); +void exposeTridiagonalization(); +void exposeComplexEigenSolver(); +void exposeComplexSchur(); void exposeLLTSolver(); void exposeLDLTSolver(); +void exposeFullPivLUSolver(); +void exposePartialPivLUSolver(); void exposeQRSolvers(); -void exposeMINRESSolver(); void exposeSimplicialLLTSolver(); void exposeSimplicialLDLTSolver(); +void exposeSparseLUSolver(); +void exposeSparseQRSolver(); void exposePermutationMatrix(); +void exposeBDCSVDSolver(); +void exposeJacobiSVDSolver(); + +void exposeBackwardCompatibilityAliases() { + typedef Eigen::MatrixXd MatrixXd; + MINRESSolverVisitor::expose("MINRES"); +} void exposeDecompositions() { using namespace Eigen; exposeEigenSolver(); + exposeGeneralizedEigenSolver(); exposeSelfAdjointEigenSolver(); + exposeGeneralizedSelfAdjointEigenSolver(); + exposeHessenbergDecomposition(); + exposeRealQZ(); + exposeRealSchur(); + exposeTridiagonalization(); + exposeComplexEigenSolver(); + exposeComplexSchur(); exposeLLTSolver(); exposeLDLTSolver(); + exposeFullPivLUSolver(); + exposePartialPivLUSolver(); exposeQRSolvers(); - exposeMINRESSolver(); + exposeBDCSVDSolver(); + exposeJacobiSVDSolver(); { bp::enum_("DecompositionOptions") @@ -44,6 +74,8 @@ void exposeDecompositions() { // Expose sparse decompositions exposeSimplicialLLTSolver(); exposeSimplicialLDLTSolver(); + exposeSparseLUSolver(); + exposeSparseQRSolver(); exposePermutationMatrix(); @@ -54,5 +86,7 @@ void exposeDecompositions() { #ifdef EIGENPY_WITH_ACCELERATE_SUPPORT exposeAccelerate(); #endif + + exposeBackwardCompatibilityAliases(); } } // namespace eigenpy diff --git a/src/decompositions/fullpivlu-solver.cpp b/src/decompositions/fullpivlu-solver.cpp new file mode 100644 index 000000000..44fcaa07d --- /dev/null +++ b/src/decompositions/fullpivlu-solver.cpp @@ -0,0 +1,12 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/FullPivLU.hpp" + +namespace eigenpy { +void exposeFullPivLUSolver() { + using namespace Eigen; + FullPivLUSolverVisitor::expose("FullPivLU"); +} +} // namespace eigenpy diff --git a/src/decompositions/generalized-eigen-solver.cpp b/src/decompositions/generalized-eigen-solver.cpp new file mode 100644 index 000000000..424393a70 --- /dev/null +++ b/src/decompositions/generalized-eigen-solver.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/GeneralizedEigenSolver.hpp" + +namespace eigenpy { +void exposeGeneralizedEigenSolver() { + using namespace Eigen; + GeneralizedEigenSolverVisitor::expose("GeneralizedEigenSolver"); +} +} // namespace eigenpy diff --git a/src/decompositions/generalized-self-adjoint-eigen-solver.cpp b/src/decompositions/generalized-self-adjoint-eigen-solver.cpp new file mode 100644 index 000000000..9dd254604 --- /dev/null +++ b/src/decompositions/generalized-self-adjoint-eigen-solver.cpp @@ -0,0 +1,14 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp" + +namespace eigenpy { +void exposeGeneralizedSelfAdjointEigenSolver() { + using namespace Eigen; + GeneralizedSelfAdjointEigenSolverVisitor::expose( + "GeneralizedSelfAdjointEigenSolver"); +} +} // namespace eigenpy diff --git a/src/decompositions/hessenberg-decomposition.cpp b/src/decompositions/hessenberg-decomposition.cpp new file mode 100644 index 000000000..6daebfc83 --- /dev/null +++ b/src/decompositions/hessenberg-decomposition.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/HessenbergDecomposition.hpp" + +namespace eigenpy { +void exposeHessenbergDecomposition() { + using namespace Eigen; + HessenbergDecompositionVisitor::expose("HessenbergDecomposition"); +} +} // namespace eigenpy diff --git a/src/decompositions/jacobisvd-solver.cpp b/src/decompositions/jacobisvd-solver.cpp new file mode 100644 index 000000000..a59538ca9 --- /dev/null +++ b/src/decompositions/jacobisvd-solver.cpp @@ -0,0 +1,29 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/JacobiSVD.hpp" + +namespace eigenpy { +void exposeJacobiSVDSolver() { + using namespace Eigen; + using Eigen::JacobiSVD; + + using Eigen::ColPivHouseholderQRPreconditioner; + using Eigen::FullPivHouseholderQRPreconditioner; + using Eigen::HouseholderQRPreconditioner; + using Eigen::NoQRPreconditioner; + + using ColPivHhJacobiSVD = + JacobiSVD; + using FullPivHhJacobiSVD = + JacobiSVD; + using HhJacobiSVD = JacobiSVD; + using NoPrecondJacobiSVD = JacobiSVD; + + JacobiSVDVisitor::expose("ColPivHhJacobiSVD"); + JacobiSVDVisitor::expose("FullPivHhJacobiSVD"); + JacobiSVDVisitor::expose("HhJacobiSVD"); + JacobiSVDVisitor::expose("NoPrecondJacobiSVD"); +} +} // namespace eigenpy diff --git a/src/decompositions/partialpivlu-solver.cpp b/src/decompositions/partialpivlu-solver.cpp new file mode 100644 index 000000000..c6ee7f9be --- /dev/null +++ b/src/decompositions/partialpivlu-solver.cpp @@ -0,0 +1,12 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/PartialPivLU.hpp" + +namespace eigenpy { +void exposePartialPivLUSolver() { + using namespace Eigen; + PartialPivLUSolverVisitor::expose("PartialPivLU"); +} +} // namespace eigenpy diff --git a/src/decompositions/real-qz.cpp b/src/decompositions/real-qz.cpp new file mode 100644 index 000000000..0ed7561e8 --- /dev/null +++ b/src/decompositions/real-qz.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/RealQZ.hpp" + +namespace eigenpy { +void exposeRealQZ() { + using namespace Eigen; + RealQZVisitor::expose("RealQZ"); +} +} // namespace eigenpy diff --git a/src/decompositions/real-schur.cpp b/src/decompositions/real-schur.cpp new file mode 100644 index 000000000..dbfe0bd7b --- /dev/null +++ b/src/decompositions/real-schur.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/RealSchur.hpp" + +namespace eigenpy { +void exposeRealSchur() { + using namespace Eigen; + RealSchurVisitor::expose("RealSchur"); +} +} // namespace eigenpy diff --git a/src/decompositions/seigen-solver.cpp b/src/decompositions/seigen-solver.cpp deleted file mode 100644 index 8fdfbf4af..000000000 --- a/src/decompositions/seigen-solver.cpp +++ /dev/null @@ -1,13 +0,0 @@ - -/* - * Copyright 2024 INRIA - */ - -#include "eigenpy/decompositions/EigenSolver.hpp" - -namespace eigenpy { -void exposeEigenSolver() { - using namespace Eigen; - EigenSolverVisitor::expose("EigenSolver"); -} -} // namespace eigenpy diff --git a/src/decompositions/simplicial-ldlt-solver.cpp b/src/decompositions/simplicial-ldlt-solver.cpp index c9a3e7b98..1d090afeb 100644 --- a/src/decompositions/simplicial-ldlt-solver.cpp +++ b/src/decompositions/simplicial-ldlt-solver.cpp @@ -2,7 +2,7 @@ * Copyright 2024 INRIA */ -#include "eigenpy/decompositions/sparse/LDLT.hpp" +#include "eigenpy/decompositions/sparse/SimplicialLDLT.hpp" namespace eigenpy { void exposeSimplicialLDLTSolver() { diff --git a/src/decompositions/simplicial-llt-solver.cpp b/src/decompositions/simplicial-llt-solver.cpp index f46e36f25..70f069970 100644 --- a/src/decompositions/simplicial-llt-solver.cpp +++ b/src/decompositions/simplicial-llt-solver.cpp @@ -2,7 +2,7 @@ * Copyright 2024 INRIA */ -#include "eigenpy/decompositions/sparse/LLT.hpp" +#include "eigenpy/decompositions/sparse/SimplicialLLT.hpp" namespace eigenpy { void exposeSimplicialLLTSolver() { diff --git a/src/decompositions/sparse-lu-solver.cpp b/src/decompositions/sparse-lu-solver.cpp new file mode 100644 index 000000000..763a75592 --- /dev/null +++ b/src/decompositions/sparse-lu-solver.cpp @@ -0,0 +1,27 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/sparse/SparseLU.hpp" + +namespace eigenpy { +void exposeSparseLUSolver() { + using namespace Eigen; + + typedef SparseMatrix ColMajorSparseMatrix; + typedef typename ColMajorSparseMatrix::StorageIndex StorageIndex; + typedef COLAMDOrdering Ordering; + typedef SparseLU SparseLUType; + + typedef typename SparseLUType::Scalar Scalar; + typedef typename SparseLUType::SCMatrix SCMatrix; + typedef Eigen::MappedSparseMatrix + MappedSparseMatrix; + + SparseLUMatrixLReturnTypeVisitor::expose( + ("SparseLUMatrixLReturnType")); + SparseLUMatrixUReturnTypeVisitor::expose( + ("SparseLUMatrixUReturnType")); + SparseLUVisitor::expose("SparseLU"); +} +} // namespace eigenpy diff --git a/src/decompositions/sparse-qr-solver.cpp b/src/decompositions/sparse-qr-solver.cpp new file mode 100644 index 000000000..0bac14197 --- /dev/null +++ b/src/decompositions/sparse-qr-solver.cpp @@ -0,0 +1,22 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/sparse/SparseQR.hpp" + +namespace eigenpy { +void exposeSparseQRSolver() { + using namespace Eigen; + + typedef SparseMatrix ColMajorSparseMatrix; + typedef typename ColMajorSparseMatrix::StorageIndex StorageIndex; + typedef COLAMDOrdering Ordering; + typedef SparseQR SparseQRType; + + SparseQRMatrixQTransposeReturnTypeVisitor::expose( + "SparseQRMatrixQTransposeReturnType"); + SparseQRMatrixQReturnTypeVisitor::expose( + "SparseQRMatrixQReturnType"); + SparseQRVisitor::expose("SparseQR"); +} +} // namespace eigenpy diff --git a/src/decompositions/tridiagonalization.cpp b/src/decompositions/tridiagonalization.cpp new file mode 100644 index 000000000..52febadd2 --- /dev/null +++ b/src/decompositions/tridiagonalization.cpp @@ -0,0 +1,13 @@ + +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/decompositions/Tridiagonalization.hpp" + +namespace eigenpy { +void exposeTridiagonalization() { + using namespace Eigen; + TridiagonalizationVisitor::expose("Tridiagonalization"); +} +} // namespace eigenpy diff --git a/src/solvers/bicgstab.cpp b/src/solvers/bicgstab.cpp new file mode 100644 index 000000000..f2fb0876f --- /dev/null +++ b/src/solvers/bicgstab.cpp @@ -0,0 +1,17 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/solvers/BiCGSTAB.hpp" + +namespace eigenpy { +void exposeBiCGSTAB() { + using namespace Eigen; + using Eigen::BiCGSTAB; + using Eigen::IdentityPreconditioner; + using IdentityBiCGSTAB = BiCGSTAB; + + BiCGSTABVisitor>::expose("BiCGSTAB"); + BiCGSTABVisitor::expose("IdentityBiCGSTAB"); +} +} // namespace eigenpy diff --git a/src/solvers/conjugate-gradient.cpp b/src/solvers/conjugate-gradient.cpp new file mode 100644 index 000000000..a64cce257 --- /dev/null +++ b/src/solvers/conjugate-gradient.cpp @@ -0,0 +1,21 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/solvers/ConjugateGradient.hpp" + +namespace eigenpy { +void exposeConjugateGradient() { + using namespace Eigen; + using Eigen::ConjugateGradient; + using Eigen::IdentityPreconditioner; + using Eigen::Lower; + using IdentityConjugateGradient = + ConjugateGradient; + + ConjugateGradientVisitor>::expose( + "ConjugateGradient"); + ConjugateGradientVisitor::expose( + "IdentityConjugateGradient"); +} +} // namespace eigenpy diff --git a/src/solvers/incomplete-cholesky.cpp b/src/solvers/incomplete-cholesky.cpp new file mode 100644 index 000000000..b6b2274de --- /dev/null +++ b/src/solvers/incomplete-cholesky.cpp @@ -0,0 +1,13 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/solvers/IncompleteCholesky.hpp" + +namespace eigenpy { +void exposeIncompleteCholesky() { + using namespace Eigen; + typedef SparseMatrix ColMajorSparseMatrix; + IncompleteCholeskyVisitor::expose("IncompleteCholesky"); +} +} // namespace eigenpy diff --git a/src/solvers/incomplete-lut.cpp b/src/solvers/incomplete-lut.cpp new file mode 100644 index 000000000..640fe6178 --- /dev/null +++ b/src/solvers/incomplete-lut.cpp @@ -0,0 +1,12 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/solvers/IncompleteLUT.hpp" + +namespace eigenpy { +void exposeIncompleteLUT() { + typedef Eigen::SparseMatrix ColMajorSparseMatrix; + IncompleteLUTVisitor::expose("IncompleteLUT"); +} +} // namespace eigenpy diff --git a/src/solvers/least-squares-conjugate-gradient.cpp b/src/solvers/least-squares-conjugate-gradient.cpp new file mode 100644 index 000000000..5fefbf67a --- /dev/null +++ b/src/solvers/least-squares-conjugate-gradient.cpp @@ -0,0 +1,29 @@ +/* + * Copyright 2025 INRIA + */ + +#include "eigenpy/solvers/LeastSquaresConjugateGradient.hpp" + +namespace eigenpy { +void exposeLeastSquaresConjugateGradient() { + using namespace Eigen; + using Eigen::LeastSquaresConjugateGradient; + + using Eigen::DiagonalPreconditioner; + using Eigen::IdentityPreconditioner; + using Eigen::LeastSquareDiagonalPreconditioner; + + using IdentityLeastSquaresConjugateGradient = + LeastSquaresConjugateGradient; + using DiagonalLeastSquaresConjugateGradient = LeastSquaresConjugateGradient< + MatrixXd, DiagonalPreconditioner>; + + LeastSquaresConjugateGradientVisitor>>:: + expose("LeastSquaresConjugateGradient"); + LeastSquaresConjugateGradientVisitor:: + expose("IdentityLeastSquaresConjugateGradient"); + LeastSquaresConjugateGradientVisitor:: + expose("DiagonalLeastSquaresConjugateGradient"); +} +} // namespace eigenpy diff --git a/src/decompositions/minres-solver.cpp b/src/solvers/minres.cpp similarity index 57% rename from src/decompositions/minres-solver.cpp rename to src/solvers/minres.cpp index 44a12cea1..06f8ff63e 100644 --- a/src/decompositions/minres-solver.cpp +++ b/src/solvers/minres.cpp @@ -1,11 +1,11 @@ /* - * Copyright 2024 INRIA + * Copyright 2025 INRIA */ -#include "eigenpy/decompositions/minres.hpp" +#include "eigenpy/solvers/MINRES.hpp" namespace eigenpy { -void exposeMINRESSolver() { +void exposeMINRES() { using namespace Eigen; MINRESSolverVisitor::expose("MINRES"); } diff --git a/src/solvers/solvers.cpp b/src/solvers/solvers.cpp index 207eb6e76..e9a0d7fd2 100644 --- a/src/solvers/solvers.cpp +++ b/src/solvers/solvers.cpp @@ -1,36 +1,26 @@ /* - * Copyright 2017-2020 CNRS INRIA + * Copyright 2017-2025 CNRS INRIA */ -#include - -#if EIGEN_VERSION_AT_LEAST(3, 2, 0) - -#include "eigenpy/solvers/ConjugateGradient.hpp" #include "eigenpy/solvers/solvers.hpp" -#if EIGEN_VERSION_AT_LEAST(3, 3, 5) -#include "eigenpy/solvers/LeastSquaresConjugateGradient.hpp" -#endif +#include "eigenpy/fwd.hpp" namespace eigenpy { -void exposeSolvers() { - using namespace Eigen; - ConjugateGradientVisitor< - ConjugateGradient>::expose(); -#if EIGEN_VERSION_AT_LEAST(3, 3, 5) - LeastSquaresConjugateGradientVisitor>>::expose(); -#endif - // Conjugate gradient with limited BFGS preconditioner - ConjugateGradientVisitor< - ConjugateGradient>:: - expose("IdentityConjugateGradient"); - // ConjugateGradientVisitor< - // ConjugateGradient - // > >::expose("LimitedBFGSConjugateGradient"); +void exposeBiCGSTAB(); +void exposeMINRES(); +void exposeConjugateGradient(); +void exposeLeastSquaresConjugateGradient(); +void exposeIncompleteCholesky(); +void exposeIncompleteLUT(); + +void exposeSolvers() { + exposeBiCGSTAB(); + exposeMINRES(); + exposeConjugateGradient(); + exposeLeastSquaresConjugateGradient(); + exposeIncompleteCholesky(); + exposeIncompleteLUT(); } } // namespace eigenpy - -#endif diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt index 98485fbed..130b3e1e2 100644 --- a/unittest/CMakeLists.txt +++ b/unittest/CMakeLists.txt @@ -117,7 +117,7 @@ add_python_lib_unit_test("py-multiple-registration" "unittest/python/test_multiple_registration.py") add_python_lib_unit_test("py-tensor" "unittest/python/test_tensor.py") -add_python_lib_unit_test("py-geometry" "unittest/python/test_geometry.py") +add_python_lib_unit_test("py-Geometry" "unittest/python/test_Geometry.py") add_python_lib_unit_test("py-complex" "unittest/python/test_complex.py") add_python_lib_unit_test("py-deprecation-policy" "unittest/python/test_deprecation_policy.py") @@ -134,21 +134,56 @@ add_python_eigenpy_lib_unit_test("py-dimensions" add_python_eigenpy_lib_unit_test("py-version" "unittest/python/test_version.py") -add_python_eigenpy_lib_unit_test("py-eigen-solver" - "unittest/python/test_eigen_solver.py") +add_python_eigenpy_lib_unit_test("py-EigenSolver" + "unittest/python/test_EigenSolver.py") add_python_eigenpy_lib_unit_test( - "py-self-adjoint-eigen-solver" - "unittest/python/test_self_adjoint_eigen_solver.py") + "py-GeneralizedEigenSolver" "unittest/python/test_GeneralizedEigenSolver.py") + +add_python_eigenpy_lib_unit_test( + "py-GeneralizedSelfAdjointEigenSolver" + "unittest/python/test_GeneralizedSelfAdjointEigenSolver.py") + +add_python_eigenpy_lib_unit_test("py-ComplexEigenSolver" + "unittest/python/test_ComplexEigenSolver.py") + +add_python_eigenpy_lib_unit_test("py-ComplexSchur" + "unittest/python/test_ComplexSchur.py") + +add_python_eigenpy_lib_unit_test( + "py-SelfAdjointEigenSolver" "unittest/python/test_SelfAdjointEigenSolver.py") + +add_python_eigenpy_lib_unit_test( + "py-HessenbergDecomposition" + "unittest/python/test_HessenbergDecomposition.py") + +add_python_eigenpy_lib_unit_test("py-RealQZ" "unittest/python/test_RealQZ.py") + +add_python_eigenpy_lib_unit_test("py-RealSchur" + "unittest/python/test_RealSchur.py") + +add_python_eigenpy_lib_unit_test("py-Tridiagonalization" + "unittest/python/test_Tridiagonalization.py") add_python_eigenpy_lib_unit_test("py-LLT" "unittest/python/test_LLT.py") add_python_eigenpy_lib_unit_test("py-LDLT" "unittest/python/test_LDLT.py") +add_python_eigenpy_lib_unit_test("py-FullPivLU" + "unittest/python/test_FullPivLU.py") + +add_python_eigenpy_lib_unit_test("py-PartialPivLU" + "unittest/python/test_PartialPivLU.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") +add_python_eigenpy_lib_unit_test("py-BDCSVD" "unittest/python/test_BDCSVD.py") + +add_python_eigenpy_lib_unit_test("py-JacobiSVD" + "unittest/python/test_JacobiSVD.py") + if(NOT WIN32) add_python_eigenpy_lib_unit_test("py-MINRES" "unittest/python/test_MINRES.py") endif(NOT WIN32) @@ -182,6 +217,14 @@ if(BUILD_TESTING_SCIPY) add_python_eigenpy_unit_test( "py-SimplicialLDLT" "unittest/python/decompositions/sparse/test_SimplicialLDLT.py") + add_python_eigenpy_unit_test("py-IncompleteCholesky" + "unittest/python/test_IncompleteCholesky.py") + add_python_eigenpy_unit_test("py-IncompleteLUT" + "unittest/python/test_IncompleteLUT.py") + add_python_eigenpy_unit_test( + "py-SparseLU" "unittest/python/decompositions/sparse/test_SparseLU.py") + add_python_eigenpy_unit_test( + "py-SparseQR" "unittest/python/decompositions/sparse/test_SparseQR.py") if(BUILD_WITH_CHOLMOD_SUPPORT) add_python_eigenpy_unit_test( diff --git a/unittest/python/decompositions/sparse/test_SparseLU.py b/unittest/python/decompositions/sparse/test_SparseLU.py new file mode 100644 index 000000000..ee0cd7f3e --- /dev/null +++ b/unittest/python/decompositions/sparse/test_SparseLU.py @@ -0,0 +1,64 @@ +import numpy as np +import scipy.sparse as spa + +import eigenpy + +dim = 100 +rng = np.random.default_rng() + +A_fac = spa.random(dim, dim, density=0.25, random_state=rng) +A = A_fac.T @ A_fac +A += spa.diags(10.0 * rng.standard_normal(dim) ** 2) +A = A.tocsc(True) +A.check_format() + +splu = eigenpy.SparseLU(A) + +assert splu.info() == eigenpy.ComputationInfo.Success + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = splu.solve(B) +assert isinstance(X_est, np.ndarray) +assert eigenpy.is_approx(X, X_est) +assert eigenpy.is_approx(A.dot(X_est), B) + +splu.analyzePattern(A) +splu.factorize(A) + +X_sparse = spa.random(dim, 10, random_state=rng) +B_sparse = A.dot(X_sparse) +B_sparse: spa.csc_matrix = B_sparse.tocsc(True) +if not B_sparse.has_sorted_indices: + B_sparse.sort_indices() + +X_est = splu.solve(B_sparse) +assert isinstance(X_est, spa.csc_matrix) +assert eigenpy.is_approx(X_est.toarray(), X_sparse.toarray()) +assert eigenpy.is_approx(A.dot(X_est.toarray()), B_sparse.toarray()) + +assert splu.nnzL() > 0 +assert splu.nnzU() > 0 + +L = splu.matrixL() +U = splu.matrixU() + +assert L.rows() == dim +assert L.cols() == dim +assert U.rows() == dim +assert U.cols() == dim + +x_true = rng.random(dim) +b_true = A.dot(x_true) +P_rows_indices = splu.rowsPermutation().indices() +P_cols_indices = splu.colsPermutation().indices() + +b_permuted = b_true[P_rows_indices] +z = b_permuted.copy() +L.solveInPlace(z) +y = z.copy() +U.solveInPlace(y) +x_reconstructed = np.zeros(dim) +x_reconstructed[P_cols_indices] = y + +assert eigenpy.is_approx(x_reconstructed, x_true, prec=1e-6) diff --git a/unittest/python/decompositions/sparse/test_SparseQR.py b/unittest/python/decompositions/sparse/test_SparseQR.py new file mode 100644 index 000000000..f6b6d070c --- /dev/null +++ b/unittest/python/decompositions/sparse/test_SparseQR.py @@ -0,0 +1,95 @@ +import numpy as np +import scipy.sparse as spa + +import eigenpy + +dim = 100 +rng = np.random.default_rng() + +A_fac = spa.random(dim, dim, density=0.25, random_state=rng) +A = A_fac.T @ A_fac +A += spa.diags(10.0 * rng.standard_normal(dim) ** 2) +A = A.tocsc(True) +A.check_format() + +spqr = eigenpy.SparseQR(A) + +assert spqr.info() == eigenpy.ComputationInfo.Success + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = spqr.solve(B) +assert isinstance(X_est, np.ndarray) +assert eigenpy.is_approx(X, X_est) +assert eigenpy.is_approx(A.dot(X_est), B) + +spqr.analyzePattern(A) +spqr.factorize(A) + +X_sparse = spa.random(dim, 10, random_state=rng) +B_sparse = A.dot(X_sparse) +B_sparse: spa.csc_matrix = B_sparse.tocsc(True) +if not B_sparse.has_sorted_indices: + B_sparse.sort_indices() + +X_est = spqr.solve(B_sparse) +assert isinstance(X_est, spa.csc_matrix) +assert eigenpy.is_approx(X_est.toarray(), X_sparse.toarray()) +assert eigenpy.is_approx(A.dot(X_est.toarray()), B_sparse.toarray()) + +Q = spqr.matrixQ() +R = spqr.matrixR() +P = spqr.colsPermutation() + +assert spqr.matrixQ().rows() == dim +assert spqr.matrixQ().cols() == dim +assert R.shape[0] == dim +assert R.shape[1] == dim +assert P.indices().size == dim + +test_vec = rng.random(dim) +test_matrix = rng.random((dim, 20)) + +Qv = Q @ test_vec +QM = Q @ test_matrix +Qt = Q.transpose() +QtV = Qt @ test_vec +QtM = Qt @ test_matrix + +assert Qv.shape == (dim,) +assert QM.shape == (dim, 20) +assert QtV.shape == (dim,) +assert QtM.shape == (dim, 20) + +Qa_real_mat = Q.adjoint() +QaV = Qa_real_mat @ test_vec +assert eigenpy.is_approx(QtV, QaV) + +A_dense = A.toarray() +P_indices = np.array([P.indices()[i] for i in range(dim)]) +A_permuted = A_dense[:, P_indices] + +QtAP = Qt @ A_permuted +R_dense = spqr.matrixR().toarray() +assert eigenpy.is_approx(QtAP, R_dense) + +Q_sparse = Q.toSparse() +R_sparse = R + +assert Q_sparse.shape == (dim, dim) + +QtQ_sparse = Q_sparse.T @ Q_sparse +QQt_sparse = Q_sparse @ Q_sparse.T +I_sparse = spa.identity(dim, format="csc") + +assert eigenpy.is_approx(QtQ_sparse.toarray(), I_sparse.toarray()) +assert eigenpy.is_approx(QQt_sparse.toarray(), I_sparse.toarray()) + +Q_sparse_test_vec = Q_sparse @ test_vec +assert eigenpy.is_approx(Qv, Q_sparse_test_vec) + +Q_sparse_test_matrix = Q_sparse @ test_matrix +assert eigenpy.is_approx(QM, Q_sparse_test_matrix) + +QR_sparse = Q_sparse @ R_sparse.toarray() +assert eigenpy.is_approx(QR_sparse, A_permuted) diff --git a/unittest/python/test_BDCSVD.py b/unittest/python/test_BDCSVD.py new file mode 100644 index 000000000..72c8e51c8 --- /dev/null +++ b/unittest/python/test_BDCSVD.py @@ -0,0 +1,138 @@ +import numpy as np + +import eigenpy + +_options = [ + 0, + eigenpy.DecompositionOptions.ComputeThinU, + eigenpy.DecompositionOptions.ComputeThinV, + eigenpy.DecompositionOptions.ComputeFullU, + eigenpy.DecompositionOptions.ComputeFullV, + eigenpy.DecompositionOptions.ComputeThinU + | eigenpy.DecompositionOptions.ComputeThinV, + eigenpy.DecompositionOptions.ComputeFullU + | eigenpy.DecompositionOptions.ComputeFullV, + eigenpy.DecompositionOptions.ComputeThinU + | eigenpy.DecompositionOptions.ComputeFullV, + eigenpy.DecompositionOptions.ComputeFullU + | eigenpy.DecompositionOptions.ComputeThinV, +] + + +def test_bdcsvd(options): + dim = 100 + rng = np.random.default_rng() + A = rng.random((dim, dim)) + A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim)) + + bdcsvd = eigenpy.BDCSVD(A, options) + assert bdcsvd.info() == eigenpy.ComputationInfo.Success + + if options & ( + eigenpy.DecompositionOptions.ComputeThinU + | eigenpy.DecompositionOptions.ComputeFullU + ) and options & ( + eigenpy.DecompositionOptions.ComputeThinV + | eigenpy.DecompositionOptions.ComputeFullV + ): + X = rng.random((dim, 20)) + B = A @ X + X_est = bdcsvd.solve(B) + assert eigenpy.is_approx(X, X_est) + assert eigenpy.is_approx(A @ X_est, B) + + x = rng.random(dim) + b = A @ x + x_est = bdcsvd.solve(b) + assert eigenpy.is_approx(x, x_est) + assert eigenpy.is_approx(A @ x_est, b) + + rows = bdcsvd.rows() + cols = bdcsvd.cols() + assert cols == dim + assert rows == dim + + _bdcsvd_compute = bdcsvd.compute(A) + _bdcsvd_compute_options = bdcsvd.compute(A, options) + + rank = bdcsvd.rank() + singularvalues = bdcsvd.singularValues() + nonzerosingularvalues = bdcsvd.nonzeroSingularValues() + assert rank == nonzerosingularvalues + assert len(singularvalues) == dim + assert all( + singularvalues[i] >= singularvalues[i + 1] + for i in range(len(singularvalues) - 1) + ) + + compute_u = bdcsvd.computeU() + compute_v = bdcsvd.computeV() + expected_compute_u = bool( + options + & ( + eigenpy.DecompositionOptions.ComputeThinU + | eigenpy.DecompositionOptions.ComputeFullU + ) + ) + expected_compute_v = bool( + options + & ( + eigenpy.DecompositionOptions.ComputeThinV + | eigenpy.DecompositionOptions.ComputeFullV + ) + ) + assert compute_u == expected_compute_u + assert compute_v == expected_compute_v + + if compute_u: + matrixU = bdcsvd.matrixU() + if options & eigenpy.DecompositionOptions.ComputeFullU: + assert matrixU.shape == (dim, dim) + elif options & eigenpy.DecompositionOptions.ComputeThinU: + assert matrixU.shape == (dim, dim) + assert eigenpy.is_approx(matrixU.T @ matrixU, np.eye(matrixU.shape[1])) + + if compute_v: + matrixV = bdcsvd.matrixV() + if options & eigenpy.DecompositionOptions.ComputeFullV: + assert matrixV.shape == (dim, dim) + elif options & eigenpy.DecompositionOptions.ComputeThinV: + assert matrixV.shape == (dim, dim) + assert eigenpy.is_approx(matrixV.T @ matrixV, np.eye(matrixV.shape[1])) + + if compute_u and compute_v: + U = bdcsvd.matrixU() + V = bdcsvd.matrixV() + S = bdcsvd.singularValues() + S_matrix = np.diag(S) + A_reconstructed = U @ S_matrix @ V.T + assert eigenpy.is_approx(A, A_reconstructed) + + bdcsvd.setSwitchSize(5) + bdcsvd.setSwitchSize(16) + bdcsvd.setSwitchSize(32) + + bdcsvd.setThreshold() + _default_threshold = bdcsvd.threshold() + bdcsvd.setThreshold(1e-8) + assert bdcsvd.threshold() == 1e-8 + + decomp1 = eigenpy.BDCSVD() + decomp2 = eigenpy.BDCSVD() + id1 = decomp1.id() + id2 = decomp2.id() + assert id1 != id2 + assert id1 == decomp1.id() + assert id2 == decomp2.id() + + decomp3 = eigenpy.BDCSVD(dim, dim, options) + decomp4 = eigenpy.BDCSVD(dim, dim, options) + id3 = decomp3.id() + id4 = decomp4.id() + assert id3 != id4 + assert id3 == decomp3.id() + assert id4 == decomp4.id() + + +for opt in _options: + test_bdcsvd(opt) diff --git a/unittest/python/test_ComplexEigenSolver.py b/unittest/python/test_ComplexEigenSolver.py new file mode 100644 index 000000000..8157a7221 --- /dev/null +++ b/unittest/python/test_ComplexEigenSolver.py @@ -0,0 +1,28 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +es = eigenpy.ComplexEigenSolver(A) +assert es.info() == eigenpy.ComputationInfo.Success + +V = es.eigenvectors() +D = es.eigenvalues() +assert V.shape == (dim, dim) +assert D.shape == (dim,) + +AV = A @ V +VD = V @ np.diag(D) +assert eigenpy.is_approx(AV.real, VD.real) +assert eigenpy.is_approx(AV.imag, VD.imag) + +ces5 = eigenpy.ComplexEigenSolver(A) +ces6 = eigenpy.ComplexEigenSolver(A) +id5 = ces5.id() +id6 = ces6.id() +assert id5 != id6 +assert id5 == ces5.id() +assert id6 == ces6.id() diff --git a/unittest/python/test_ComplexSchur.py b/unittest/python/test_ComplexSchur.py new file mode 100644 index 000000000..f777dc930 --- /dev/null +++ b/unittest/python/test_ComplexSchur.py @@ -0,0 +1,48 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +cs = eigenpy.ComplexSchur(A) +assert cs.info() == eigenpy.ComputationInfo.Success + +U = cs.matrixU() +T = cs.matrixT() + +A_complex = A.astype(complex) +assert eigenpy.is_approx(A_complex, U @ T @ U.conj().T) +assert eigenpy.is_approx(U @ U.conj().T, np.eye(dim)) + +for row in range(1, dim): + for col in range(row): + assert abs(T[row, col]) < 1e-12 + +A_triangular = np.triu(A) +cs_triangular = eigenpy.ComplexSchur(dim) +cs_triangular.setMaxIterations(1) +result_triangular = cs_triangular.compute(A_triangular) +assert result_triangular.info() == eigenpy.ComputationInfo.Success + +T_triangular = cs_triangular.matrixT() +U_triangular = cs_triangular.matrixU() + +A_triangular_complex = A_triangular.astype(complex) +assert eigenpy.is_approx(T_triangular, A_triangular_complex) +assert eigenpy.is_approx(U_triangular, np.eye(dim, dtype=complex)) + +hess = eigenpy.HessenbergDecomposition(A) +H = hess.matrixH() +Q_hess = hess.matrixQ() + +cs_from_hess = eigenpy.ComplexSchur(dim) +result_from_hess = cs_from_hess.computeFromHessenberg(H, Q_hess, True) +assert result_from_hess.info() == eigenpy.ComputationInfo.Success + +T_from_hess = cs_from_hess.matrixT() +U_from_hess = cs_from_hess.matrixU() + +A_complex = A.astype(complex) +assert eigenpy.is_approx(A_complex, U_from_hess @ T_from_hess @ U_from_hess.conj().T) diff --git a/unittest/python/test_eigen_solver.py b/unittest/python/test_EigenSolver.py similarity index 100% rename from unittest/python/test_eigen_solver.py rename to unittest/python/test_EigenSolver.py diff --git a/unittest/python/test_FullPivLU.py b/unittest/python/test_FullPivLU.py new file mode 100644 index 000000000..cdbbc4ad4 --- /dev/null +++ b/unittest/python/test_FullPivLU.py @@ -0,0 +1,109 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim)) +fullpivlu = eigenpy.FullPivLU(A) + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = fullpivlu.solve(B) +assert eigenpy.is_approx(X, X_est) +assert eigenpy.is_approx(A.dot(X_est), B) + +x = rng.random(dim) +b = A.dot(x) +x_est = fullpivlu.solve(b) +assert eigenpy.is_approx(x, x_est) +assert eigenpy.is_approx(A.dot(x_est), b) + +rows = fullpivlu.rows() +cols = fullpivlu.cols() +assert cols == dim +assert rows == dim + +fullpivlu_compute = fullpivlu.compute(A) +A_reconstructed = fullpivlu.reconstructedMatrix() +assert eigenpy.is_approx(A_reconstructed, A) + +nonzeropivots = fullpivlu.nonzeroPivots() +maxpivot = fullpivlu.maxPivot() +assert nonzeropivots == dim +assert maxpivot > 0 + +LU = fullpivlu.matrixLU() +P_perm = fullpivlu.permutationP() +Q_perm = fullpivlu.permutationQ() +P = P_perm.toDenseMatrix() +Q = Q_perm.toDenseMatrix() + +U = np.triu(LU) +L = np.eye(dim) + np.tril(LU, -1) +assert eigenpy.is_approx(P @ A @ Q, L @ U) + +rank = fullpivlu.rank() +dimkernel = fullpivlu.dimensionOfKernel() +injective = fullpivlu.isInjective() +surjective = fullpivlu.isSurjective() +invertible = fullpivlu.isInvertible() +assert rank == dim +assert dimkernel == 0 +assert injective +assert surjective +assert invertible + +kernel = fullpivlu.kernel() +assert eigenpy.is_approx(A @ kernel, np.zeros((dim, 1)), 1e-10) +image = fullpivlu.image(A) +assert image.shape[1] == rank + +inverse = fullpivlu.inverse() +assert eigenpy.is_approx(A @ inverse, np.eye(dim)) +assert eigenpy.is_approx(inverse @ A, np.eye(dim)) + +rcond = fullpivlu.rcond() +determinant = fullpivlu.determinant() +det_numpy = np.linalg.det(A) +assert rcond > 0 +assert abs(determinant - det_numpy) / abs(det_numpy) < 1e-10 + +fullpivlu.setThreshold() +default_threshold = fullpivlu.threshold() +fullpivlu.setThreshold(1e-8) +assert fullpivlu.threshold() == 1e-8 + +P_inv = P_perm.inverse().toDenseMatrix() +Q_inv = Q_perm.inverse().toDenseMatrix() +assert eigenpy.is_approx(P @ P_inv, np.eye(dim)) +assert eigenpy.is_approx(Q @ Q_inv, np.eye(dim)) +assert eigenpy.is_approx(P_inv @ P, np.eye(dim)) +assert eigenpy.is_approx(Q_inv @ Q, np.eye(dim)) + +rows_rect = 4 +cols_rect = 6 +A_rect = rng.random((rows_rect, cols_rect)) +fullpivlu_rect = eigenpy.FullPivLU(A_rect) +assert fullpivlu_rect.rows() == rows_rect +assert fullpivlu_rect.cols() == cols_rect +rank_rect = fullpivlu_rect.rank() +assert rank_rect <= min(rows_rect, cols_rect) +assert fullpivlu_rect.dimensionOfKernel() == cols_rect - rank_rect + +decomp1 = eigenpy.FullPivLU() +decomp2 = eigenpy.FullPivLU() +id1 = decomp1.id() +id2 = decomp2.id() +assert id1 != id2 +assert id1 == decomp1.id() +assert id2 == decomp2.id() + +decomp3 = eigenpy.FullPivLU(dim, dim) +decomp4 = eigenpy.FullPivLU(dim, dim) +id3 = decomp3.id() +id4 = decomp4.id() +assert id3 != id4 +assert id3 == decomp3.id() +assert id4 == decomp4.id() diff --git a/unittest/python/test_GeneralizedEigenSolver.py b/unittest/python/test_GeneralizedEigenSolver.py new file mode 100644 index 000000000..66748e5fb --- /dev/null +++ b/unittest/python/test_GeneralizedEigenSolver.py @@ -0,0 +1,41 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +B = rng.random((dim, dim)) +B = (B + B.T) * 0.5 + np.diag(10.0 + rng.random(dim)) + +ges_matrices = eigenpy.GeneralizedEigenSolver(A, B) +assert ges_matrices.info() == eigenpy.ComputationInfo.Success + +alphas = ges_matrices.alphas() +betas = ges_matrices.betas() +eigenvectors = ges_matrices.eigenvectors() +eigenvalues = ges_matrices.eigenvalues() + +for k in range(dim): + v = eigenvectors[:, k] + lambda_k = eigenvalues[k] + + Av = A @ v + lambda_Bv = lambda_k * (B @ v) + assert eigenpy.is_approx(Av.real, lambda_Bv.real, 1e-6) + assert eigenpy.is_approx(Av.imag, lambda_Bv.imag, 1e-6) + +for k in range(dim): + v = eigenvectors[:, k] + alpha = alphas[k] + beta = betas[k] + + alpha_Bv = alpha * (B @ v) + beta_Av = beta * (A @ v) + assert eigenpy.is_approx(alpha_Bv.real, beta_Av.real, 1e-6) + assert eigenpy.is_approx(alpha_Bv.imag, beta_Av.imag, 1e-6) + +for k in range(dim): + if abs(betas[k]) > 1e-12: + expected_eigenvalue = alphas[k] / betas[k] + assert abs(eigenvalues[k] - expected_eigenvalue) < 1e-12 diff --git a/unittest/python/test_GeneralizedSelfAdjointEigenSolver.py b/unittest/python/test_GeneralizedSelfAdjointEigenSolver.py new file mode 100644 index 000000000..78ea15564 --- /dev/null +++ b/unittest/python/test_GeneralizedSelfAdjointEigenSolver.py @@ -0,0 +1,79 @@ +import numpy as np + +import eigenpy + +_options = [ + eigenpy.DecompositionOptions.ComputeEigenvectors + | eigenpy.DecompositionOptions.Ax_lBx, + eigenpy.DecompositionOptions.EigenvaluesOnly | eigenpy.DecompositionOptions.Ax_lBx, + eigenpy.DecompositionOptions.ComputeEigenvectors + | eigenpy.DecompositionOptions.ABx_lx, + eigenpy.DecompositionOptions.EigenvaluesOnly | eigenpy.DecompositionOptions.ABx_lx, + eigenpy.DecompositionOptions.ComputeEigenvectors + | eigenpy.DecompositionOptions.BAx_lx, + eigenpy.DecompositionOptions.EigenvaluesOnly | eigenpy.DecompositionOptions.BAx_lx, +] + + +def test_generalized_selfadjoint_eigensolver(options): + dim = 100 + rng = np.random.default_rng() + A = rng.random((dim, dim)) + A = (A + A.T) * 0.5 + B = rng.random((dim, dim)) + B = B @ B.T + 0.1 * np.eye(dim) + + gsaes = eigenpy.GeneralizedSelfAdjointEigenSolver(A, B, options) + assert gsaes.info() == eigenpy.ComputationInfo.Success + + D = gsaes.eigenvalues() + assert D.shape == (dim,) + assert all(abs(D[i].imag) < 1e-12 for i in range(dim)) + assert all(D[i] <= D[i + 1] + 1e-12 for i in range(dim - 1)) + + compute_eigenvectors = bool( + options & eigenpy.DecompositionOptions.ComputeEigenvectors + ) + + if compute_eigenvectors: + V = gsaes.eigenvectors() + assert V.shape == (dim, dim) + + if options & eigenpy.DecompositionOptions.Ax_lBx: + for i in range(dim): + v = V[:, i] + lam = D[i] + Av = A @ v + lam_Bv = lam * (B @ v) + assert eigenpy.is_approx(Av, lam_Bv, 1e-6) + + VT_B_V = V.T @ B @ V + assert eigenpy.is_approx(VT_B_V, np.eye(dim), 1e-6) + + elif options & eigenpy.DecompositionOptions.ABx_lx: + AB = A @ B + for i in range(dim): + v = V[:, i] + lam = D[i] + ABv = AB @ v + lam_v = lam * v + assert eigenpy.is_approx(ABv, lam_v, 1e-6) + + elif options & eigenpy.DecompositionOptions.BAx_lx: + BA = B @ A + for i in range(dim): + v = V[:, i] + lam = D[i] + BAv = BA @ v + lam_v = lam * v + assert eigenpy.is_approx(BAv, lam_v, 1e-6) + + _gsaes_compute = gsaes.compute(A, B) + _gsaes_compute_options = gsaes.compute(A, B, options) + + rank = len([d for d in D if abs(d) > 1e-12]) + assert rank <= dim + + +for opt in _options: + test_generalized_selfadjoint_eigensolver(opt) diff --git a/unittest/python/test_geometry.py b/unittest/python/test_Geometry.py similarity index 100% rename from unittest/python/test_geometry.py rename to unittest/python/test_Geometry.py diff --git a/unittest/python/test_HessenbergDecomposition.py b/unittest/python/test_HessenbergDecomposition.py new file mode 100644 index 000000000..095ac7105 --- /dev/null +++ b/unittest/python/test_HessenbergDecomposition.py @@ -0,0 +1,71 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +hess = eigenpy.HessenbergDecomposition(A) + +Q = hess.matrixQ() +H = hess.matrixH() + +if np.iscomplexobj(A): + A_reconstructed = Q @ H @ Q.conj().T +else: + A_reconstructed = Q @ H @ Q.T +assert eigenpy.is_approx(A, A_reconstructed) + +for row in range(2, dim): + for col in range(row - 1): + assert abs(H[row, col]) < 1e-12 + +if np.iscomplexobj(Q): + QQ_conj = Q @ Q.conj().T +else: + QQ_conj = Q @ Q.T +assert eigenpy.is_approx(QQ_conj, np.eye(dim)) + +A_test = rng.random((dim, dim)) +hess1 = eigenpy.HessenbergDecomposition(dim) +hess1.compute(A_test) +hess2 = eigenpy.HessenbergDecomposition(A_test) + +H1 = hess1.matrixH() +H2 = hess2.matrixH() +Q1 = hess1.matrixQ() +Q2 = hess2.matrixQ() + +assert eigenpy.is_approx(H1, H2) +assert eigenpy.is_approx(Q1, Q2) + +hCoeffs = hess.householderCoefficients() +packed = hess.packedMatrix() + +assert hCoeffs.shape == (dim - 1,) +assert packed.shape == (dim, dim) + +for i in range(dim): + for j in range(i - 1, dim): + if j >= 0: + assert abs(H[i, j] - packed[i, j]) < 1e-12 + +hess_default = eigenpy.HessenbergDecomposition(dim) +hess_matrix = eigenpy.HessenbergDecomposition(A) + +hess1_id = eigenpy.HessenbergDecomposition(dim) +hess2_id = eigenpy.HessenbergDecomposition(dim) +id1 = hess1_id.id() +id2 = hess2_id.id() +assert id1 != id2 +assert id1 == hess1_id.id() +assert id2 == hess2_id.id() + +hess3_id = eigenpy.HessenbergDecomposition(A) +hess4_id = eigenpy.HessenbergDecomposition(A) +id3 = hess3_id.id() +id4 = hess4_id.id() +assert id3 != id4 +assert id3 == hess3_id.id() +assert id4 == hess4_id.id() diff --git a/unittest/python/test_IncompleteCholesky.py b/unittest/python/test_IncompleteCholesky.py new file mode 100644 index 000000000..b634f8b90 --- /dev/null +++ b/unittest/python/test_IncompleteCholesky.py @@ -0,0 +1,71 @@ +import numpy as np +from scipy.sparse import csc_matrix + +import eigenpy + +dim = 100 +rng = np.random.default_rng() + +A = rng.random((dim, dim)) +A = (A + A.T) * 0.5 + np.diag(5.0 + rng.random(dim)) +A = csc_matrix(A) + +ichol = eigenpy.solvers.IncompleteCholesky(A) +assert ichol.info() == eigenpy.ComputationInfo.Success +assert ichol.rows() == dim +assert ichol.cols() == dim + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = ichol.solve(B) +assert isinstance(X_est, np.ndarray) +residual = np.linalg.norm(B - A.dot(X_est)) / np.linalg.norm(B) +assert residual < 0.1 + +x = rng.random(dim) +b = A.dot(x) +x_est = ichol.solve(b) +assert isinstance(x_est, np.ndarray) +residual = np.linalg.norm(b - A.dot(x_est)) / np.linalg.norm(b) +assert residual < 0.1 + +X_sparse = csc_matrix(rng.random((dim, 10))) +B_sparse = A.dot(X_sparse).tocsc() +if not B_sparse.has_sorted_indices: + B_sparse.sort_indices() +X_est_sparse = ichol.solve(B_sparse) +assert isinstance(X_est_sparse, csc_matrix) + +ichol.analyzePattern(A) +ichol.factorize(A) +ichol.compute(A) +assert ichol.info() == eigenpy.ComputationInfo.Success + +L = ichol.matrixL() +S_diag = ichol.scalingS() +perm = ichol.permutationP() +P = perm.toDenseMatrix() + +assert isinstance(L, csc_matrix) +assert isinstance(S_diag, np.ndarray) +assert L.shape == (dim, dim) +assert S_diag.shape == (dim,) + +L_dense = L.toarray() +upper_part = np.triu(L_dense, k=1) +assert np.allclose(upper_part, 0, atol=1e-12) + +assert np.all(S_diag > 0) + +S = csc_matrix((S_diag, (range(dim), range(dim))), shape=(dim, dim)) + +PA = P @ A +PAP = PA @ P.T +SPAP = S @ PAP +SPAPS = SPAP @ S + +LLT = L @ L.T + +diff = SPAPS - LLT +relative_error = np.linalg.norm(diff.data) / np.linalg.norm(SPAPS.data) +assert relative_error < 0.5 diff --git a/unittest/python/test_IncompleteLUT.py b/unittest/python/test_IncompleteLUT.py new file mode 100644 index 000000000..0ff820ef7 --- /dev/null +++ b/unittest/python/test_IncompleteLUT.py @@ -0,0 +1,50 @@ +import numpy as np +from scipy.sparse import csc_matrix + +import eigenpy + +dim = 100 +rng = np.random.default_rng() + +A = rng.random((dim, dim)) +A = (A + A.T) * 0.5 + np.diag(5.0 + rng.random(dim)) +A = csc_matrix(A) + +ilut = eigenpy.solvers.IncompleteLUT(A) +assert ilut.info() == eigenpy.ComputationInfo.Success +assert ilut.rows() == dim +assert ilut.cols() == dim + +X = rng.random((dim, 100)) +B = A.dot(X) +X_est = ilut.solve(B) +assert isinstance(X_est, np.ndarray) +residual = np.linalg.norm(B - A.dot(X_est)) / np.linalg.norm(B) +assert residual < 0.1 + +x = rng.random(dim) +b = A.dot(x) +x_est = ilut.solve(b) +assert isinstance(x_est, np.ndarray) +residual = np.linalg.norm(b - A.dot(x_est)) / np.linalg.norm(b) +assert residual < 0.1 + +X_sparse = csc_matrix(rng.random((dim, 10))) +B_sparse = A.dot(X_sparse).tocsc() +if not B_sparse.has_sorted_indices: + B_sparse.sort_indices() +X_est_sparse = ilut.solve(B_sparse) +assert isinstance(X_est_sparse, csc_matrix) + +ilut.analyzePattern(A) +ilut.factorize(A) +assert ilut.info() == eigenpy.ComputationInfo.Success + +ilut_params = eigenpy.solvers.IncompleteLUT(A, 1e-4, 15) +assert ilut_params.info() == eigenpy.ComputationInfo.Success + +ilut_set = eigenpy.solvers.IncompleteLUT() +ilut_set.setDroptol(1e-3) +ilut_set.setFillfactor(20) +ilut_set.compute(A) +assert ilut_set.info() == eigenpy.ComputationInfo.Success diff --git a/unittest/python/test_JacobiSVD.py b/unittest/python/test_JacobiSVD.py new file mode 100644 index 000000000..a9c975c90 --- /dev/null +++ b/unittest/python/test_JacobiSVD.py @@ -0,0 +1,114 @@ +import numpy as np + +import eigenpy + +THIN_U = eigenpy.DecompositionOptions.ComputeThinU +THIN_V = eigenpy.DecompositionOptions.ComputeThinV +FULL_U = eigenpy.DecompositionOptions.ComputeFullU +FULL_V = eigenpy.DecompositionOptions.ComputeFullV + +_options = [ + 0, + # THIN_U, + # THIN_V, + # FULL_U, + # FULL_V, + THIN_U | THIN_V, + FULL_U | FULL_V, + # THIN_U | FULL_V, + # FULL_U | THIN_V, +] + +_classes = [ + eigenpy.ColPivHhJacobiSVD, + # eigenpy.FullPivHhJacobiSVD, + # eigenpy.HhJacobiSVD, + # eigenpy.NoPrecondJacobiSVD, +] + +# Rationale: Tets only few cases to gain computation time +# User can test all of them by uncommenting the corresponding lines + + +def is_valid_combination(cls, opt): + if cls == eigenpy.FullPivHhJacobiSVD: + has_thin_u = bool(opt & THIN_U) + has_thin_v = bool(opt & THIN_V) + + if has_thin_u or has_thin_v: + return False + + return True + + +def test_jacobi(cls, opt): + dim = 100 + rng = np.random.default_rng() + A = rng.random((dim, dim)) + A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim)) + + jacobisvd = cls(A, opt) + assert jacobisvd.info() == eigenpy.ComputationInfo.Success + + has_u = opt & (THIN_U | FULL_U) + has_v = opt & (THIN_V | FULL_V) + + if has_u and has_v: + X = rng.random((dim, 20)) + B = A @ X + X_est = jacobisvd.solve(B) + assert eigenpy.is_approx(X, X_est) + assert eigenpy.is_approx(A @ X_est, B) + + x = rng.random(dim) + b = A @ x + x_est = jacobisvd.solve(b) + assert eigenpy.is_approx(x, x_est) + assert eigenpy.is_approx(A @ x_est, b) + + assert jacobisvd.rows() == dim + assert jacobisvd.cols() == dim + + _jacobisvd_compute = jacobisvd.compute(A) + _jacobisvd_compute_options = jacobisvd.compute(A, opt) + + rank = jacobisvd.rank() + singularvalues = jacobisvd.singularValues() + nonzerosingularvalues = jacobisvd.nonzeroSingularValues() + assert rank == nonzerosingularvalues + assert len(singularvalues) == dim + assert all( + singularvalues[i] >= singularvalues[i + 1] + for i in range(len(singularvalues) - 1) + ) + + compute_u = jacobisvd.computeU() + compute_v = jacobisvd.computeV() + expected_compute_u = bool(has_u) + expected_compute_v = bool(has_v) + assert compute_u == expected_compute_u + assert compute_v == expected_compute_v + + if compute_u: + matrixU = jacobisvd.matrixU() + assert matrixU.shape == (dim, dim) + assert eigenpy.is_approx(matrixU.T @ matrixU, np.eye(matrixU.shape[1])) + + if compute_v: + matrixV = jacobisvd.matrixV() + assert matrixV.shape == (dim, dim) + assert eigenpy.is_approx(matrixV.T @ matrixV, np.eye(matrixV.shape[1])) + + if compute_u and compute_v: + U = jacobisvd.matrixU() + V = jacobisvd.matrixV() + S = jacobisvd.singularValues() + S_matrix = np.diag(S) + A_reconstructed = U @ S_matrix @ V.T + assert eigenpy.is_approx(A, A_reconstructed) + + +for cls in _classes: + for opt in _options: + if is_valid_combination(cls, opt): + test_jacobi(cls, opt) diff --git a/unittest/python/test_MINRES.py b/unittest/python/test_MINRES.py index b630ad4bc..c830a4f77 100644 --- a/unittest/python/test_MINRES.py +++ b/unittest/python/test_MINRES.py @@ -6,11 +6,16 @@ rng = np.random.default_rng() A = np.eye(dim) -minres = eigenpy.MINRES(A) +minres = eigenpy.solvers.MINRES(A) X = rng.random((dim, 20)) B = A.dot(X) X_est = minres.solve(B) -print("A.dot(X_est):", A.dot(X_est)) -print("B:", B) +assert eigenpy.is_approx(A.dot(X_est), B, 1e-6) + +minres_back = eigenpy.MINRES(A) + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = minres_back.solve(B) assert eigenpy.is_approx(A.dot(X_est), B, 1e-6) diff --git a/unittest/python/test_PartialPivLU.py b/unittest/python/test_PartialPivLU.py new file mode 100644 index 000000000..487a515f5 --- /dev/null +++ b/unittest/python/test_PartialPivLU.py @@ -0,0 +1,76 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim)) +partialpivlu = eigenpy.PartialPivLU(A) + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = partialpivlu.solve(B) +assert eigenpy.is_approx(X, X_est) +assert eigenpy.is_approx(A.dot(X_est), B) + +x = rng.random(dim) +b = A.dot(x) +x_est = partialpivlu.solve(b) +assert eigenpy.is_approx(x, x_est) +assert eigenpy.is_approx(A.dot(x_est), b) + +rows = partialpivlu.rows() +cols = partialpivlu.cols() +assert cols == dim +assert rows == dim + +partialpivlu_compute = partialpivlu.compute(A) +A_reconstructed = partialpivlu.reconstructedMatrix() +assert eigenpy.is_approx(A_reconstructed, A) + +LU = partialpivlu.matrixLU() +P_perm = partialpivlu.permutationP() +P = P_perm.toDenseMatrix() + +U = np.triu(LU) +L = np.eye(dim) + np.tril(LU, -1) +assert eigenpy.is_approx(P @ A, L @ U) + +inverse = partialpivlu.inverse() +assert eigenpy.is_approx(A @ inverse, np.eye(dim)) +assert eigenpy.is_approx(inverse @ A, np.eye(dim)) + +rcond = partialpivlu.rcond() +determinant = partialpivlu.determinant() +det_numpy = np.linalg.det(A) +assert rcond > 0 +assert abs(determinant - det_numpy) / abs(det_numpy) < 1e-10 + +P_inv = P_perm.inverse().toDenseMatrix() +assert eigenpy.is_approx(P @ P_inv, np.eye(dim)) +assert eigenpy.is_approx(P_inv @ P, np.eye(dim)) + +decomp1 = eigenpy.PartialPivLU() +decomp2 = eigenpy.PartialPivLU() +id1 = decomp1.id() +id2 = decomp2.id() +assert id1 != id2 +assert id1 == decomp1.id() +assert id2 == decomp2.id() + +decomp3 = eigenpy.PartialPivLU(dim) +decomp4 = eigenpy.PartialPivLU(dim) +id3 = decomp3.id() +id4 = decomp4.id() +assert id3 != id4 +assert id3 == decomp3.id() +assert id4 == decomp4.id() + +decomp5 = eigenpy.PartialPivLU(A) +decomp6 = eigenpy.PartialPivLU(A) +id5 = decomp5.id() +id6 = decomp6.id() +assert id5 != id6 +assert id5 == decomp5.id() +assert id6 == decomp6.id() diff --git a/unittest/python/test_RealQZ.py b/unittest/python/test_RealQZ.py new file mode 100644 index 000000000..df6262d43 --- /dev/null +++ b/unittest/python/test_RealQZ.py @@ -0,0 +1,38 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +B = rng.random((dim, dim)) + +realqz = eigenpy.RealQZ(A, B) +assert realqz.info() == eigenpy.ComputationInfo.Success + +Q = realqz.matrixQ() +S = realqz.matrixS() +Z = realqz.matrixZ() +T = realqz.matrixT() + +assert eigenpy.is_approx(A, Q @ S @ Z) +assert eigenpy.is_approx(B, Q @ T @ Z) + +assert eigenpy.is_approx(Q @ Q.T, np.eye(dim)) +assert eigenpy.is_approx(Z @ Z.T, np.eye(dim)) + +for i in range(dim): + for j in range(i): + assert abs(T[i, j]) < 1e-12 + +for i in range(dim): + for j in range(i - 1): + assert abs(S[i, j]) < 1e-12 + +realqz3_id = eigenpy.RealQZ(A, B) +realqz4_id = eigenpy.RealQZ(A, B) +id3 = realqz3_id.id() +id4 = realqz4_id.id() +assert id3 != id4 +assert id3 == realqz3_id.id() +assert id4 == realqz4_id.id() diff --git a/unittest/python/test_RealSchur.py b/unittest/python/test_RealSchur.py new file mode 100644 index 000000000..069d41394 --- /dev/null +++ b/unittest/python/test_RealSchur.py @@ -0,0 +1,49 @@ +import numpy as np + +import eigenpy + + +def verify_is_quasi_triangular(T): + size = T.shape[0] + + for row in range(2, size): + for col in range(row - 1): + assert abs(T[row, col]) < 1e-12 + + for row in range(1, size): + if abs(T[row, row - 1]) > 1e-12: + if row < size - 1: + assert abs(T[row + 1, row]) < 1e-12 + + tr = T[row - 1, row - 1] + T[row, row] + det = T[row - 1, row - 1] * T[row, row] - T[row - 1, row] * T[row, row - 1] + assert 4 * det > tr * tr + + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +rs = eigenpy.RealSchur(A) +assert rs.info() == eigenpy.ComputationInfo.Success + +U = rs.matrixU() +T = rs.matrixT() + +assert eigenpy.is_approx(A, U @ T @ U.T) +assert eigenpy.is_approx(U @ U.T, np.eye(dim)) + +verify_is_quasi_triangular(T) + +hess = eigenpy.HessenbergDecomposition(A) +H = hess.matrixH() +Q_hess = hess.matrixQ() + +rs_from_hess = eigenpy.RealSchur(dim) +result_from_hess = rs_from_hess.computeFromHessenberg(H, Q_hess, True) +assert result_from_hess.info() == eigenpy.ComputationInfo.Success + +T_from_hess = rs_from_hess.matrixT() +U_from_hess = rs_from_hess.matrixU() + +assert eigenpy.is_approx(A, U_from_hess @ T_from_hess @ U_from_hess.T) diff --git a/unittest/python/test_self_adjoint_eigen_solver.py b/unittest/python/test_SelfAdjointEigenSolver.py similarity index 100% rename from unittest/python/test_self_adjoint_eigen_solver.py rename to unittest/python/test_SelfAdjointEigenSolver.py diff --git a/unittest/python/test_Tridiagonalization.py b/unittest/python/test_Tridiagonalization.py new file mode 100644 index 000000000..1270c778d --- /dev/null +++ b/unittest/python/test_Tridiagonalization.py @@ -0,0 +1,104 @@ +import numpy as np + +import eigenpy + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +A = (A + A.T) * 0.5 + +tri = eigenpy.Tridiagonalization(A) + +Q = tri.matrixQ() +T = tri.matrixT() + +assert eigenpy.is_approx(A, Q @ T @ Q.T) +assert eigenpy.is_approx(Q @ Q.T, np.eye(dim)) + +for i in range(dim): + for j in range(dim): + if abs(i - j) > 1: + assert abs(T[i, j]) < 1e-12 + +assert eigenpy.is_approx(T, T.T) + +diag = tri.diagonal() +sub_diag = tri.subDiagonal() + +for i in range(dim): + assert abs(diag[i] - T[i, i]) < 1e-12 + +for i in range(dim - 1): + assert abs(sub_diag[i] - T[i + 1, i]) < 1e-12 + +A_test = rng.random((dim, dim)) +A_test = (A_test + A_test.T) * 0.5 + +tri1 = eigenpy.Tridiagonalization(dim) +tri1.compute(A_test) +tri2 = eigenpy.Tridiagonalization(A_test) + +Q1 = tri1.matrixQ() +T1 = tri1.matrixT() +Q2 = tri2.matrixQ() +T2 = tri2.matrixT() + +assert eigenpy.is_approx(Q1, Q2) +assert eigenpy.is_approx(T1, T2) + +h_coeffs = tri.householderCoefficients() +packed = tri.packedMatrix() + +assert h_coeffs.shape == (dim - 1,) +assert packed.shape == (dim, dim) + +for i in range(dim): + for j in range(i + 1, dim): + assert abs(packed[i, j] - A[i, j]) < 1e-12 + +for i in range(dim): + assert abs(packed[i, i] - T[i, i]) < 1e-12 + if i < dim - 1: + assert abs(packed[i + 1, i] - T[i + 1, i]) < 1e-12 + +A_diag = np.diag(rng.random(dim)) +tri_diag = eigenpy.Tridiagonalization(A_diag) +Q_diag = tri_diag.matrixQ() +T_diag = tri_diag.matrixT() + +assert eigenpy.is_approx(A_diag, Q_diag @ T_diag @ Q_diag.T) +for i in range(dim): + for j in range(dim): + if i != j: + assert abs(T_diag[i, j]) < 1e-10 + +A_tridiag = np.zeros((dim, dim)) +for i in range(dim): + A_tridiag[i, i] = rng.random() + if i < dim - 1: + val = rng.random() + A_tridiag[i, i + 1] = val + A_tridiag[i + 1, i] = val + +tri_tridiag = eigenpy.Tridiagonalization(A_tridiag) +Q_tridiag = tri_tridiag.matrixQ() +T_tridiag = tri_tridiag.matrixT() + +assert eigenpy.is_approx(A_tridiag, Q_tridiag @ T_tridiag @ Q_tridiag.T) + + +tri1_id = eigenpy.Tridiagonalization(dim) +tri2_id = eigenpy.Tridiagonalization(dim) +id1 = tri1_id.id() +id2 = tri2_id.id() +assert id1 != id2 +assert id1 == tri1_id.id() +assert id2 == tri2_id.id() + +tri3_id = eigenpy.Tridiagonalization(A) +tri4_id = eigenpy.Tridiagonalization(A) +id3 = tri3_id.id() +id4 = tri4_id.id() +assert id3 != id4 +assert id3 == tri3_id.id() +assert id4 == tri4_id.id()