diff --git a/CMakeLists.txt b/CMakeLists.txt index ed255dd..3af56f2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -131,7 +131,7 @@ target_include_directories( ) target_link_libraries(nanoeigenpy_headers INTERFACE Eigen3::Eigen) -set(${PROJECT_NAME}_SOURCES src/module.cpp src/solvers.cpp) +set(${PROJECT_NAME}_SOURCES src/module.cpp) nanobind_add_module(nanoeigenpy NB_STATIC NB_SUPPRESS_WARNINGS ${nanoeigenpy_SOURCES} ${nanoeigenpy_HEADERS}) target_link_libraries(nanoeigenpy PRIVATE nanoeigenpy_headers) diff --git a/cmake b/cmake index 9982101..2ede15d 160000 --- a/cmake +++ b/cmake @@ -1 +1 @@ -Subproject commit 9982101b22135509d684e2228dac747e7afa3bcf +Subproject commit 2ede15d1cb9d66401ba96788444ad64c44ffccd8 diff --git a/include/nanoeigenpy/constants.hpp b/include/nanoeigenpy/constants.hpp index 0cdb1d2..06753c2 100644 --- a/include/nanoeigenpy/constants.hpp +++ b/include/nanoeigenpy/constants.hpp @@ -30,6 +30,14 @@ inline void exposeConstants(nb::module_ m) { ._c(Ax_lBx) ._c(ABx_lx) ._c(BAx_lx); +#undef _c + using Eigen::TransformTraits; +#define _c(name) value(#name, TransformTraits::name) + nb::enum_(m, "TransformTraits") + ._c(Isometry) + ._c(Affine) + ._c(AffineCompact) + ._c(Projective); #undef _c } } // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions.hpp b/include/nanoeigenpy/decompositions.hpp index f16b582..85a63cb 100644 --- a/include/nanoeigenpy/decompositions.hpp +++ b/include/nanoeigenpy/decompositions.hpp @@ -10,10 +10,24 @@ #include "nanoeigenpy/decompositions/complete-orthogonal-decomposition.hpp" #include "nanoeigenpy/decompositions/eigen-solver.hpp" #include "nanoeigenpy/decompositions/self-adjoint-eigen-solver.hpp" +#include "nanoeigenpy/decompositions/generalized-self-adjoint-eigen-solver.hpp" +#include "nanoeigenpy/decompositions/complex-eigen-solver.hpp" +#include "nanoeigenpy/decompositions/complex-schur.hpp" +#include "nanoeigenpy/decompositions/generalized-eigen-solver.hpp" +#include "nanoeigenpy/decompositions/hessenberg-decomposition.hpp" +#include "nanoeigenpy/decompositions/real-qz.hpp" +#include "nanoeigenpy/decompositions/real-schur.hpp" +#include "nanoeigenpy/decompositions/tridiagonalization.hpp" #include "nanoeigenpy/decompositions/permutation-matrix.hpp" +#include "nanoeigenpy/decompositions/full-piv-lu.hpp" +#include "nanoeigenpy/decompositions/partial-piv-lu.hpp" +#include "nanoeigenpy/decompositions/bdcsvd.hpp" +#include "nanoeigenpy/decompositions/jacobi-svd.hpp" #include "nanoeigenpy/decompositions/sparse/simplicial-llt.hpp" #include "nanoeigenpy/decompositions/sparse/simplicial-ldlt.hpp" +#include "nanoeigenpy/decompositions/sparse/sparse-lu.hpp" +#include "nanoeigenpy/decompositions/sparse/sparse-qr.hpp" #ifdef NANOEIGENPY_HAS_CHOLMOD #include "nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-llt.hpp" diff --git a/include/nanoeigenpy/decompositions/bdcsvd.hpp b/include/nanoeigenpy/decompositions/bdcsvd.hpp new file mode 100644 index 0000000..5559104 --- /dev/null +++ b/include/nanoeigenpy/decompositions/bdcsvd.hpp @@ -0,0 +1,89 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "nanoeigenpy/decompositions/svd-base.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +MatrixOrVector solve(const Eigen::BDCSVD &c, + const MatrixOrVector &vec) { + return c.solve(vec); +} + +template +void exposeBDCSVD(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::BDCSVD; + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + nb::class_( + m, name, + "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.") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), + "rows"_a, "cols"_a, "computationOptions"_a = 0, + "Default constructor with memory preallocation.") + + .def(nb::init(), "matrix"_a, + "computationOptions"_a = 0, + "Constructs a SVD factorization from a given matrix.") + + .def(SVDBaseVisitor()) + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes the SVD of given matrix.", + nb::rv_policy::reference) + .def( + "compute", + [](Solver &c, const MatrixType &matrix, unsigned int) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "computationOptions"_a, + "Computes the SVD of given matrix.", nb::rv_policy::reference) + + .def("setSwitchSize", &Solver::setSwitchSize, "s"_a) + + .def( + "solve", + [](const Solver &c, const VectorType &b) -> VectorType { + return solve(c, b); + }, + "b"_a, + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def( + "solve", + [](const Solver &c, const MatrixType &B) -> MatrixType { + return solve(c, B); + }, + "B"_a, + "Returns the solution X of A X = B using the current " + "decomposition of A where B is a right hand side matrix.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/col-piv-householder-qr.hpp b/include/nanoeigenpy/decompositions/col-piv-householder-qr.hpp index 7c95e8e..40bf7f7 100644 --- a/include/nanoeigenpy/decompositions/col-piv-householder-qr.hpp +++ b/include/nanoeigenpy/decompositions/col-piv-householder-qr.hpp @@ -3,11 +3,11 @@ #pragma once #include "nanoeigenpy/fwd.hpp" -#include "nanoeigenpy/eigen-base.hpp" #include namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template MatrixOrVector solve(const Eigen::ColPivHouseholderQR &c, @@ -21,7 +21,7 @@ MatrixType inverse(const Eigen::ColPivHouseholderQR &c) { } template -void exposeColPivHouseholderQRSolver(nb::module_ m, const char *name) { +void exposeColPivHouseholderQR(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Solver = Eigen::ColPivHouseholderQR; using Scalar = typename MatrixType::Scalar; @@ -48,12 +48,11 @@ void exposeColPivHouseholderQRSolver(nb::module_ m, const char *name) { "The default constructor is useful in cases in which the " "user intends to perform decompositions via " "HouseholderQR.compute(matrix).") - .def(nb::init(), nb::arg("rows"), - nb::arg("cols"), + .def(nb::init(), "rows"_a, "cols"_a, "Default constructor with memory preallocation.\n" "Like the default constructor but with preallocation of the " "internal data according to the specified problem size. ") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a QR factorization from a given matrix.\n" "This constructor computes the QR factorization of the matrix " "matrix by calling the method compute().") @@ -128,10 +127,10 @@ void exposeColPivHouseholderQRSolver(nb::module_ m, const char *name) { .def( "setThreshold", - [](Solver &c, RealScalar const &threshold) { + [](Solver &c, const RealScalar &threshold) { return c.setThreshold(threshold); }, - nb::arg("threshold"), + "threshold"_a, "Allows to prescribe a threshold to be used by certain methods, " "such as rank(), who need to determine when pivots are to be " "considered nonzero. This is not used for the QR decomposition " @@ -147,6 +146,12 @@ void exposeColPivHouseholderQRSolver(nb::module_ m, const char *name) { "is strictly greater than |pivot| ⩽ threshold×|maxpivot| where " "maxpivot is the biggest pivot.", nb::rv_policy::reference) + .def( + "setThreshold", + [](Solver &c) { return c.setThreshold(Eigen::Default); }, + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + nb::rv_policy::reference) .def("threshold", &Solver::threshold, "Returns the threshold that will be used by certain methods such " "as rank().") @@ -161,31 +166,31 @@ void exposeColPivHouseholderQRSolver(nb::module_ m, const char *name) { .def( "compute", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the QR factorization of given matrix.", + "matrix"_a, "Computes the QR factorization of given matrix.", nb::rv_policy::reference) .def( - "inverse", [](Solver const &c) -> MatrixType { return inverse(c); }, + "inverse", [](const Solver &c) -> MatrixType { return inverse(c); }, "Returns the inverse of the matrix associated with the QR " "decomposition.") .def( "solve", - [](Solver const &c, VectorType const &b) -> VectorType { + [](const Solver &c, const VectorType &b) -> VectorType { return solve(c, b); }, - nb::arg("b"), + "b"_a, "Returns the solution x of A x = B using the current " "decomposition of A where b is a right hand side vector.") .def( "solve", - [](Solver const &c, MatrixType const &B) -> MatrixType { + [](const Solver &c, const MatrixType &B) -> MatrixType { return solve(c, B); }, - nb::arg("B"), + "B"_a, "Returns the solution X of A X = B using the current " "decomposition of A where B is a right hand side matrix.") diff --git a/include/nanoeigenpy/decompositions/complete-orthogonal-decomposition.hpp b/include/nanoeigenpy/decompositions/complete-orthogonal-decomposition.hpp index 0ed54eb..edc1954 100644 --- a/include/nanoeigenpy/decompositions/complete-orthogonal-decomposition.hpp +++ b/include/nanoeigenpy/decompositions/complete-orthogonal-decomposition.hpp @@ -3,11 +3,11 @@ #pragma once #include "nanoeigenpy/fwd.hpp" -#include "nanoeigenpy/eigen-base.hpp" #include namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template MatrixOrVector solve( @@ -23,8 +23,7 @@ MatrixType pseudoInverse( } template -void exposeCompleteOrthogonalDecompositionSolver(nb::module_ m, - const char *name) { +void exposeCompleteOrthogonalDecomposition(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Solver = Eigen::CompleteOrthogonalDecomposition; using Scalar = typename MatrixType::Scalar; @@ -50,12 +49,11 @@ void exposeCompleteOrthogonalDecompositionSolver(nb::module_ m, "The default constructor is useful in cases in which the " "user intends to perform decompositions via " "HouseholderQR.compute(matrix).") - .def(nb::init(), nb::arg("rows"), - nb::arg("cols"), + .def(nb::init(), "rows"_a, "cols"_a, "Default constructor with memory preallocation.\n" "Like the default constructor but with preallocation of the " "internal data according to the specified problem size. ") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a QR factorization from a given matrix.\n" "This constructor computes the QR factorization of the matrix " "matrix by calling the method compute().") @@ -137,10 +135,10 @@ void exposeCompleteOrthogonalDecompositionSolver(nb::module_ m, .def( "setThreshold", - [](Solver &c, RealScalar const &threshold) { + [](Solver &c, const RealScalar &threshold) { return c.setThreshold(threshold); }, - nb::arg("threshold"), + "threshold"_a, "Allows to prescribe a threshold to be used by certain methods, " "such as rank(), who need to determine when pivots are to be " "considered nonzero. This is not used for the complete " @@ -158,6 +156,12 @@ void exposeCompleteOrthogonalDecompositionSolver(nb::module_ m, "is strictly greater than |pivot| ⩽ threshold×|maxpivot| where " "maxpivot is the biggest pivot.", nb::rv_policy::reference) + .def( + "setThreshold", + [](Solver &c) { return c.setThreshold(Eigen::Default); }, + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + nb::rv_policy::reference) .def("threshold", &Solver::threshold, "Returns the threshold that will be used by certain methods such " "as rank().") @@ -174,32 +178,32 @@ void exposeCompleteOrthogonalDecompositionSolver(nb::module_ m, .def( "compute", - [](Solver &c, MatrixType const &matrix) { return c.compute(matrix); }, - nb::arg("matrix"), + [](Solver &c, const MatrixType &matrix) { return c.compute(matrix); }, + "matrix"_a, "Computes the complete orthogonal factorization of given matrix.", nb::rv_policy::reference) .def( "pseudoInverse", - [](Solver const &c) -> MatrixType { return pseudoInverse(c); }, + [](const Solver &c) -> MatrixType { return pseudoInverse(c); }, "Returns the pseudo-inverse of the matrix associated with the " "complete orthogonal " "decomposition.") .def( "solve", - [](Solver const &c, VectorType const &b) -> VectorType { + [](const Solver &c, const VectorType &b) -> VectorType { return solve(c, b); }, - nb::arg("b"), + "b"_a, "Returns the solution x of A x = B using the current " "decomposition of A where b is a right hand side vector.") .def( "solve", - [](Solver const &c, MatrixType const &B) -> MatrixType { + [](const Solver &c, const MatrixType &B) -> MatrixType { return solve(c, B); }, - nb::arg("B"), + "B"_a, "Returns the solution X of A X = B using the current " "decomposition of A where B is a right hand side matrix.") diff --git a/include/nanoeigenpy/decompositions/complex-eigen-solver.hpp b/include/nanoeigenpy/decompositions/complex-eigen-solver.hpp new file mode 100644 index 0000000..0dd5343 --- /dev/null +++ b/include/nanoeigenpy/decompositions/complex-eigen-solver.hpp @@ -0,0 +1,64 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeComplexEigenSolver(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::ComplexEigenSolver; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Complex Eigen Solver") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "computeEigenvectors"_a = true, + "Computes eigendecomposition of given matrix") + + .def("eigenvalues", &Solver::eigenvalues, + "Returns the eigenvalues of given matrix.", + nb::rv_policy::reference_internal) + .def("eigenvectors", &Solver::eigenvectors, + "Returns the eigenvectors of given matrix.", + nb::rv_policy::reference_internal) + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes the eigendecomposition of given matrix.", + nb::rv_policy::reference) + .def( + "compute", + [](Solver &c, const MatrixType &matrix, bool computeEigenvectors) + -> Solver & { return c.compute(matrix, computeEigenvectors); }, + "matrix"_a, "computeEigenvectors"_a, + "Computes the eigendecomposition of given matrix.", + nb::rv_policy::reference) + + .def("info", &Solver::info, + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("setMaxIterations", &Solver::setMaxIterations, + "Sets the maximum number of iterations allowed.", + nb::rv_policy::reference) + .def("getMaxIterations", &Solver::getMaxIterations, + "Returns the maximum number of iterations.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/complex-schur.hpp b/include/nanoeigenpy/decompositions/complex-schur.hpp new file mode 100644 index 0000000..1ae6505 --- /dev/null +++ b/include/nanoeigenpy/decompositions/complex-schur.hpp @@ -0,0 +1,74 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeComplexSchur(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::ComplexSchur; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Complex Schur decomposition") + + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "computeU"_a = true, + "Constructor; computes Schur decomposition of given matrix.") + + .def("matrixU", &Solver::matrixU, + "Returns the unitary matrix in the Schur decomposition.", + nb::rv_policy::reference_internal) + .def("matrixT", &Solver::matrixT, + "Returns the triangular matrix in the Schur decomposition. ", + nb::rv_policy::reference_internal) + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes Schur decomposition of given matrix. ", + nb::rv_policy::reference) + .def( + "compute", + [](Solver &c, const MatrixType &matrix, bool computeU) -> Solver & { + return c.compute(matrix, computeU); + }, + "matrix"_a, "computeU"_a, + "Computes Schur decomposition of given matrix. ", + nb::rv_policy::reference) + + .def( + "computeFromHessenberg", + [](Solver &c, const MatrixType &matrixH, const MatrixType &matrixQ, + bool computeU) -> Solver & { + return c.computeFromHessenberg(matrixH, matrixQ, computeU); + }, + "matrixH"_a, "matrixQ"_a, "computeU"_a, + "Compute Schur decomposition from a given Hessenberg matrix. ", + nb::rv_policy::reference) + + .def("info", &Solver::info, + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("getMaxIterations", &Solver::getMaxIterations, + "Returns the maximum number of iterations.") + .def("setMaxIterations", &Solver::setMaxIterations, + "Sets the maximum number of iterations allowed.", + nb::rv_policy::reference) + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/eigen-solver.hpp b/include/nanoeigenpy/decompositions/eigen-solver.hpp index 7c6f880..96630dd 100644 --- a/include/nanoeigenpy/decompositions/eigen-solver.hpp +++ b/include/nanoeigenpy/decompositions/eigen-solver.hpp @@ -3,12 +3,11 @@ #pragma once #include "nanoeigenpy/fwd.hpp" -#include "nanoeigenpy/eigen-base.hpp" -#include #include namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template void exposeEigenSolver(nb::module_ m, const char *name) { @@ -20,10 +19,10 @@ void exposeEigenSolver(nb::module_ m, const char *name) { } nb::class_(m, name, "Eigen solver.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("size"), + .def(nb::init(), "size"_a, "Default constructor with memory preallocation.") - .def(nb::init(), nb::arg("matrix"), - nb::arg("compute_eigen_vectors") = true, + .def(nb::init(), "matrix"_a, + "compute_eigen_vectors"_a = true, "Computes eigendecomposition of given matrix") .def("eigenvalues", &Solver::eigenvalues, @@ -34,16 +33,16 @@ void exposeEigenSolver(nb::module_ m, const char *name) { .def( "compute", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the eigendecomposition of given matrix.", + "matrix"_a, "Computes the eigendecomposition of given matrix.", nb::rv_policy::reference) .def( "compute", - [](Solver &c, MatrixType const &matrix, bool compute_eigen_vectors) + [](Solver &c, const MatrixType &matrix, bool compute_eigen_vectors) -> Solver & { return c.compute(matrix, compute_eigen_vectors); }, - nb::arg("matrix"), nb::arg("compute_eigen_vectors"), + "matrix"_a, "compute_eigen_vectors"_a, "Computes the eigendecomposition of given matrix.", nb::rv_policy::reference) diff --git a/include/nanoeigenpy/decompositions/full-piv-householder-qr.hpp b/include/nanoeigenpy/decompositions/full-piv-householder-qr.hpp index eb71a12..3a53491 100644 --- a/include/nanoeigenpy/decompositions/full-piv-householder-qr.hpp +++ b/include/nanoeigenpy/decompositions/full-piv-householder-qr.hpp @@ -3,11 +3,11 @@ #pragma once #include "nanoeigenpy/fwd.hpp" -#include "nanoeigenpy/eigen-base.hpp" #include namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template MatrixOrVector solve(const Eigen::FullPivHouseholderQR &c, @@ -21,7 +21,7 @@ MatrixType inverse(const Eigen::FullPivHouseholderQR &c) { } template -void exposeFullPivHouseholderQRSolver(nb::module_ m, const char *name) { +void exposeFullPivHouseholderQR(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Solver = Eigen::FullPivHouseholderQR; using Scalar = typename MatrixType::Scalar; @@ -51,12 +51,11 @@ void exposeFullPivHouseholderQRSolver(nb::module_ m, const char *name) { "The default constructor is useful in cases in which the " "user intends to perform decompositions via " "HouseholderQR.compute(matrix).") - .def(nb::init(), nb::arg("rows"), - nb::arg("cols"), + .def(nb::init(), "rows"_a, "cols"_a, "Default constructor with memory preallocation.\n" "Like the default constructor but with preallocation of the " "internal data according to the specified problem size. ") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a QR factorization from a given matrix.\n" "This constructor computes the QR factorization of the matrix " "matrix by calling the method compute().") @@ -126,10 +125,10 @@ void exposeFullPivHouseholderQRSolver(nb::module_ m, const char *name) { .def( "setThreshold", - [](Solver &c, RealScalar const &threshold) { + [](Solver &c, const RealScalar &threshold) { return c.setThreshold(threshold); }, - nb::arg("threshold"), + "threshold"_a, "Allows to prescribe a threshold to be used by certain methods, " "such as rank(), who need to determine when pivots are to be " "considered nonzero. This is not used for the QR decomposition " @@ -145,6 +144,12 @@ void exposeFullPivHouseholderQRSolver(nb::module_ m, const char *name) { "is strictly greater than |pivot| ⩽ threshold×|maxpivot| where " "maxpivot is the biggest pivot.", nb::rv_policy::reference) + .def( + "setThreshold", + [](Solver &c) { return c.setThreshold(Eigen::Default); }, + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + nb::rv_policy::reference) .def("threshold", &Solver::threshold, "Returns the threshold that will be used by certain methods such " "as rank().") @@ -156,31 +161,31 @@ void exposeFullPivHouseholderQRSolver(nb::module_ m, const char *name) { .def( "compute", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the QR factorization of given matrix.", + "matrix"_a, "Computes the QR factorization of given matrix.", nb::rv_policy::reference) .def( - "inverse", [](Solver const &c) -> MatrixType { return inverse(c); }, + "inverse", [](const Solver &c) -> MatrixType { return inverse(c); }, "Returns the inverse of the matrix associated with the QR " "decomposition.") .def( "solve", - [](Solver const &c, VectorType const &b) -> VectorType { + [](const Solver &c, const VectorType &b) -> VectorType { return solve(c, b); }, - nb::arg("b"), + "b"_a, "Returns the solution x of A x = B using the current " "decomposition of A where b is a right hand side vector.") .def( "solve", - [](Solver const &c, MatrixType const &B) -> MatrixType { + [](const Solver &c, const MatrixType &B) -> MatrixType { return solve(c, B); }, - nb::arg("B"), + "B"_a, "Returns the solution X of A X = B using the current " "decomposition of A where B is a right hand side matrix.") diff --git a/include/nanoeigenpy/decompositions/full-piv-lu.hpp b/include/nanoeigenpy/decompositions/full-piv-lu.hpp new file mode 100644 index 0000000..e369b6c --- /dev/null +++ b/include/nanoeigenpy/decompositions/full-piv-lu.hpp @@ -0,0 +1,207 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +MatrixOrVector solve(const Eigen::FullPivLU &c, + const MatrixOrVector &vec) { + return c.solve(vec); +} + +template +void exposeFullPivLU(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::FullPivLU; + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + using VectorType = Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, + "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.\n\n" + "The data of the LU decomposition can be directly " + "accessed through the methods " + "matrixLU(), permutationP(), permutationQ().") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "rows"_a, "cols"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "Constructs a LU factorization from a given matrix.") + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes the LU of given matrix.", + nb::rv_policy::reference) + + .def("matrixLU", &Solver::matrixLU, + "Returns the LU decomposition matrix: the upper-triangular part is " + "U, the " + "unit-lower-triangular part is L (at least for square matrices; in " + "the non-square " + "case, special care is needed, see the documentation of class " + "FullPivLU).", + nb::rv_policy::reference_internal) + + .def("nonzeroPivots", &Solver::nonzeroPivots, + "Returns the number of nonzero pivots in the LU decomposition.") + .def("maxPivot", &Solver::maxPivot, + "Returns the absolute value of the biggest pivot, i.e. the biggest" + "diagonal coefficient of U.") + + .def("permutationP", &Solver::permutationP, + "Returns the permutation matrix P in the decomposition A = P^{-1} L " + "U Q^{-1}.", + nb::rv_policy::reference_internal) + .def("permutationQ", &Solver::permutationQ, + "Returns the permutation matrix Q in the decomposition A = P^{-1} L " + "U Q^{-1}.", + nb::rv_policy::reference_internal) + + .def( + "kernel", [](Solver &c) -> MatrixType { return c.kernel(); }, + "Computes the LU of given matrix.") + .def( + "image", + [](Solver &c, const MatrixType &originalMatrix) -> MatrixType { + return c.image(originalMatrix); + }, + "Computes the LU of given matrix.") + + .def("rcond", &Solver::rcond, + "Returns an estimate of the reciprocal condition number of the " + "matrix.") + .def("determinant", &Solver::determinant, + "Returns the determinant of the underlying matrix from the " + "current factorization.") + + .def( + "setThreshold", + [](Solver &c, const RealScalar &threshold) { + return c.setThreshold(threshold); + }, + "threshold"_a, + "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.", + nb::rv_policy::reference) + .def( + "setThreshold", + [](Solver &c) { return c.setThreshold(Eigen::Default); }, + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + nb::rv_policy::reference) + .def("threshold", &Solver::threshold, + "Returns the threshold that will be used by certain methods such " + "as rank().") + + .def("rank", &Solver::rank, + "Returns the rank of the matrix associated with the LU " + "decomposition.\n" + "\n" + "Note: This method has to determine which pivots should be " + "considered nonzero. For that, it uses the threshold value that " + "you can control by calling setThreshold(threshold).") + .def("dimensionOfKernel", &Solver::dimensionOfKernel, + "Returns the dimension of the kernel of the matrix of which " + "*this is the LU decomposition.") + + .def("isInjective", &Solver::isInjective, + "Returns true if the matrix of which *this is the LU decomposition " + "represents an injective linear map, i.e. has trivial kernel; " + "false otherwise.\n\n" + "Note: This method has to determine which pivots should be " + "considered nonzero. For that, it uses the threshold value that " + "you can control by calling setThreshold(threshold).") + .def("isSurjective", &Solver::isSurjective, + "Returns true if the matrix of which *this is the LU decomposition " + "represents a surjective linear map; false otherwise.\n\n" + "Note: This method has to determine which pivots should be " + "considered nonzero. For that, it uses the threshold value that " + "you can control by calling setThreshold(threshold).") + .def("isInvertible", &Solver::isInvertible, + "Returns true if the matrix of which *this is the LU decomposition " + "is invertible.\n\n" + "Note: This method has to determine which pivots should be " + "considered nonzero. For that, it uses the threshold value that " + "you can control by calling setThreshold(threshold).") + + .def( + "inverse", [](const Solver &c) -> MatrixType { return c.inverse(); }, + "Returns the inverse of the matrix associated with the LU " + "decomposition.") + .def("reconstructedMatrix", &Solver::reconstructedMatrix, + "Returns the matrix represented by the decomposition," + "i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$." + "This function is provided for debug purposes.") + + .def("rows", &Solver::rows, "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, "Returns the number of cols of the matrix.") + + .def( + "solve", + [](const Solver &c, const VectorType &b) -> VectorType { + return solve(c, b); + }, + "b"_a, + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def( + "solve", + [](const Solver &c, const MatrixType &B) -> MatrixType { + return solve(c, B); + }, + "B"_a, + "Returns the solution X of A X = B using the current " + "decomposition of A where B is a right hand side matrix.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/generalized-eigen-solver.hpp b/include/nanoeigenpy/decompositions/generalized-eigen-solver.hpp new file mode 100644 index 0000000..e4b2c47 --- /dev/null +++ b/include/nanoeigenpy/decompositions/generalized-eigen-solver.hpp @@ -0,0 +1,70 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeGeneralizedEigenSolver(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::GeneralizedEigenSolver; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Generalized Eigen solver") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "A"_a, + "B"_a, "computeEigenvectors"_a = true, + "Constructor; computes the generalized eigendecomposition of given " + "matrix pair.") + + .def("eigenvectors", &Solver::eigenvectors, + "Returns the computed generalized eigenvectors.") + .def( + "eigenvalues", [](const Solver &c) { return c.eigenvalues().eval(); }, + "Returns the computed generalized eigenvalues.") + + .def("alphas", &Solver::alphas, + "Returns the vectors containing the alpha values.") + .def("betas", &Solver::betas, + "Returns tthe vectors containing the beta values.") + + .def( + "compute", + [](Solver &c, const MatrixType &A, const MatrixType &B) -> Solver & { + return c.compute(A, B); + }, + "A"_a, "B"_a, + "Computes generalized eigendecomposition of given matrix.", + nb::rv_policy::reference) + .def( + "compute", + [](Solver &c, const MatrixType &A, const MatrixType &B, + bool computeEigenvectors) -> Solver & { + return c.compute(A, B, computeEigenvectors); + }, + "A"_a, "B"_a, "computeEigenvectors"_a, + "Computes generalized eigendecomposition of given matrix.", + nb::rv_policy::reference) + + .def("info", &Solver::info, + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def("setMaxIterations", &Solver::setMaxIterations, + "Sets the maximum number of iterations allowed.", + nb::rv_policy::reference) + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/generalized-self-adjoint-eigen-solver.hpp b/include/nanoeigenpy/decompositions/generalized-self-adjoint-eigen-solver.hpp new file mode 100644 index 0000000..0478127 --- /dev/null +++ b/include/nanoeigenpy/decompositions/generalized-self-adjoint-eigen-solver.hpp @@ -0,0 +1,92 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeGeneralizedSelfAdjointEigenSolver(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::GeneralizedSelfAdjointEigenSolver; + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Generalized self adjoint Eigen Solver") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matA"_a, + "matB"_a, "options"_a = Eigen::ComputeEigenvectors | Eigen::Ax_lBx, + "Computes the generalized eigendecomposition of given matrix pencil") + + .def( + "compute", + [](Solver &c, const MatrixType &matA, const MatrixType &matB) + -> Solver & { return c.compute(matA, matB); }, + "matA"_a, "matB"_a, + "Computes the generalized eigendecomposition of given matrix.", + nb::rv_policy::reference) + .def( + "compute", + [](Solver &c, const MatrixType &matA, const MatrixType &matB, + int options) -> Solver & { + return c.compute(matA, matB, options); + }, + "matA"_a, "matB"_a, "options"_a, + "Computes the generalized eigendecomposition of given matrix.", + nb::rv_policy::reference) + + .def( + "eigenvalues", + [](const Solver &c) -> const VectorType & { return c.eigenvalues(); }, + "Returns the eigenvalues of given matrix.", + nb::rv_policy::reference_internal) + .def( + "eigenvectors", + [](const Solver &c) -> const MatrixType & { + return c.eigenvectors(); + }, + "Returns the eigenvectors of given matrix.", + nb::rv_policy::reference_internal) + + .def( + "computeDirect", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return static_cast(c.computeDirect(matrix)); + }, + "matrix"_a, + "Computes eigendecomposition of given matrix using a closed-form " + "algorithm.", + nb::rv_policy::reference) + .def( + "computeDirect", + [](Solver &c, const MatrixType &matrix, int options) -> Solver & { + return static_cast(c.computeDirect(matrix, options)); + }, + "matrix"_a, "options"_a, + "Computes eigendecomposition of given matrix using a closed-form " + "algorithm.", + nb::rv_policy::reference) + + .def("operatorInverseSqrt", &Solver::operatorInverseSqrt, + "Computes the inverse square root of the matrix.") + .def("operatorSqrt", &Solver::operatorSqrt, + "Computes the square root of the matrix.") + + .def("info", &Solver::info, + "NumericalIssue if the input contains INF or NaN values or " + "overflow occured. Returns Success otherwise.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/hessenberg-decomposition.hpp b/include/nanoeigenpy/decompositions/hessenberg-decomposition.hpp new file mode 100644 index 0000000..0661ed1 --- /dev/null +++ b/include/nanoeigenpy/decompositions/hessenberg-decomposition.hpp @@ -0,0 +1,53 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeHessenbergDecomposition(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::HessenbergDecomposition; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Hessenberg decomposition") + + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "Constructor; computes Hessenberg decomposition of given matrix.") + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes Hessenberg decomposition of given matrix.", + nb::rv_policy::reference) + + .def("householderCoefficients", &Solver::householderCoefficients, + "Returns the Householder coefficients.", + nb::rv_policy::reference_internal) + + .def("packedMatrix", &Solver::packedMatrix, + "Returns the internal representation of the decomposition.", + nb::rv_policy::reference_internal) + + .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(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/householder-qr.hpp b/include/nanoeigenpy/decompositions/householder-qr.hpp index 53590b5..b0cd4db 100644 --- a/include/nanoeigenpy/decompositions/householder-qr.hpp +++ b/include/nanoeigenpy/decompositions/householder-qr.hpp @@ -3,11 +3,11 @@ #pragma once #include "nanoeigenpy/fwd.hpp" -#include "nanoeigenpy/eigen-base.hpp" #include namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template MatrixOrVector solve(const Eigen::HouseholderQR &c, @@ -16,7 +16,7 @@ MatrixOrVector solve(const Eigen::HouseholderQR &c, } template -void exposeHouseholderQRSolver(nb::module_ m, const char *name) { +void exposeHouseholderQR(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Solver = Eigen::HouseholderQR; using Scalar = typename MatrixType::Scalar; @@ -48,12 +48,11 @@ void exposeHouseholderQRSolver(nb::module_ m, const char *name) { "The default constructor is useful in cases in which the " "user intends to perform decompositions via " "HouseholderQR.compute(matrix).") - .def(nb::init(), nb::arg("rows"), - nb::arg("cols"), + .def(nb::init(), "rows"_a, "cols"_a, "Default constructor with memory preallocation.\n" "Like the default constructor but with preallocation of the " "internal data according to the specified problem size. ") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a QR factorization from a given matrix.\n" "This constructor computes the QR factorization of the matrix " "matrix by calling the method compute().") @@ -89,26 +88,26 @@ void exposeHouseholderQRSolver(nb::module_ m, const char *name) { .def( "compute", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the QR factorization of given matrix.", + "matrix"_a, "Computes the QR factorization of given matrix.", nb::rv_policy::reference) .def( "solve", - [](Solver const &c, VectorType const &b) -> VectorType { + [](const Solver &c, const VectorType &b) -> VectorType { return solve(c, b); }, - nb::arg("b"), + "b"_a, "Returns the solution x of A x = B using the current " "decomposition of A where b is a right hand side vector.") .def( "solve", - [](Solver const &c, MatrixType const &B) -> MatrixType { + [](const Solver &c, const MatrixType &B) -> MatrixType { return solve(c, B); }, - nb::arg("B"), + "B"_a, "Returns the solution X of A X = B using the current " "decomposition of A where B is a right hand side matrix.") diff --git a/include/nanoeigenpy/decompositions/jacobi-svd.hpp b/include/nanoeigenpy/decompositions/jacobi-svd.hpp new file mode 100644 index 0000000..c56cae0 --- /dev/null +++ b/include/nanoeigenpy/decompositions/jacobi-svd.hpp @@ -0,0 +1,86 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "nanoeigenpy/decompositions/svd-base.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +MatrixOrVector solve(const JacobiSVD &c, const MatrixOrVector &vec) { + return c.solve(vec); +} + +template +struct JacobiSVDVisitor : nb::def_visitor> { + using MatrixType = typename JacobiSVD::MatrixType; + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Matrix; + + template + void execute(nb::class_ &cl) { + using namespace nb::literals; + cl.def(nb::init<>(), "Default constructor.") + .def(nb::init(), + "rows"_a, "cols"_a, "computationOptions"_a = 0, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "computationOptions"_a = 0, + "Constructs a SVD factorization from a given matrix.") + + .def(SVDBaseVisitor()) + + .def( + "compute", + [](JacobiSVD &c, const MatrixType &matrix) -> JacobiSVD & { + return c.compute(matrix); + }, + "matrix"_a, "Computes the SVD of given matrix.", + nb::rv_policy::reference) + .def( + "compute", + [](JacobiSVD &c, const MatrixType &matrix, + unsigned int computationOptions) -> JacobiSVD & { + return c.compute(matrix, computationOptions); + }, + "matrix"_a, "computationOptions"_a, + "Computes the SVD of given matrix.", nb::rv_policy::reference) + + .def( + "solve", + [](const JacobiSVD &c, const VectorType &b) -> VectorType { + return solve(c, b); + }, + "b"_a, + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def( + "solve", + [](const JacobiSVD &c, const MatrixType &B) -> MatrixType { + return solve(c, B); + }, + "B"_a, + "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(nb::module_ &m, const char *name) { + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name) + .def(JacobiSVDVisitor()) + .def(IdVisitor()); + } +}; + +template +void exposeJacobiSVD(nb::module_ &m, const char *name) { + JacobiSVDVisitor::expose(m, name); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/ldlt.hpp b/include/nanoeigenpy/decompositions/ldlt.hpp index 9cb8b91..f1072f4 100644 --- a/include/nanoeigenpy/decompositions/ldlt.hpp +++ b/include/nanoeigenpy/decompositions/ldlt.hpp @@ -8,6 +8,7 @@ namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template MatrixOrVector solve(const Eigen::LDLT &c, @@ -16,7 +17,7 @@ MatrixOrVector solve(const Eigen::LDLT &c, } template -void exposeLDLTSolver(nb::module_ m, const char *name) { +void exposeLDLT(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Solver = Eigen::LDLT; using Scalar = typename MatrixType::Scalar; @@ -40,9 +41,9 @@ void exposeLDLTSolver(nb::module_ m, const char *name) { "square root on D also stabilizes the computation.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("size"), + .def(nb::init(), "size"_a, "Default constructor with memory preallocation.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LLT factorization from a given matrix.") .def(EigenBaseVisitor()) @@ -53,13 +54,13 @@ void exposeLDLTSolver(nb::module_ m, const char *name) { "Returns true if the matrix is positive (semidefinite).") .def( - "matrixL", [](Solver const &c) -> MatrixType { return c.matrixL(); }, + "matrixL", [](const Solver &c) -> MatrixType { return c.matrixL(); }, "Returns the lower triangular matrix L.") .def( - "matrixU", [](Solver const &c) -> MatrixType { return c.matrixU(); }, + "matrixU", [](const Solver &c) -> MatrixType { return c.matrixU(); }, "Returns the upper triangular matrix U.") .def( - "vectorD", [](Solver const &c) -> VectorType { return c.vectorD(); }, + "vectorD", [](const Solver &c) -> VectorType { return c.vectorD(); }, "Returns the coefficients of the diagonal matrix D.") .def("matrixLDLT", &Solver::matrixLDLT, "Returns the LDLT decomposition matrix made of the lower matrix " @@ -69,7 +70,7 @@ void exposeLDLTSolver(nb::module_ m, const char *name) { .def( "transpositionsP", - [](Solver const &c) -> MatrixType { + [](const Solver &c) -> MatrixType { return c.transpositionsP() * MatrixType::Identity(c.matrixL().rows(), c.matrixL().rows()); }, @@ -77,11 +78,10 @@ void exposeLDLTSolver(nb::module_ m, const char *name) { .def( "rankUpdate", - [](Solver &c, VectorType const &w, Scalar sigma) -> Solver & { + [](Solver &c, const VectorType &w, Scalar sigma) -> Solver & { return c.rankUpdate(w, sigma); }, - "If LDL^* = A, then it becomes A + sigma * v v^*", nb::arg("w"), - nb::arg("sigma")) + "If LDL^* = A, then it becomes A + sigma * v v^*", "w"_a, "sigma"_a) .def("adjoint", &Solver::adjoint, "Returns the adjoint, that is, a reference to the decomposition " @@ -90,10 +90,10 @@ void exposeLDLTSolver(nb::module_ m, const char *name) { .def( "compute", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the LDLT of given matrix.", + "matrix"_a, "Computes the LDLT of given matrix.", nb::rv_policy::reference) .def("info", &Solver::info, "NumericalIssue if the input contains INF or NaN values or " @@ -110,18 +110,18 @@ void exposeLDLTSolver(nb::module_ m, const char *name) { .def( "solve", - [](Solver const &c, VectorType const &b) -> VectorType { + [](const Solver &c, const VectorType &b) -> VectorType { return solve(c, b); }, - nb::arg("b"), + "b"_a, "Returns the solution x of A x = b using the current " "decomposition of A.") .def( "solve", - [](Solver const &c, MatrixType const &B) -> MatrixType { + [](const Solver &c, const MatrixType &B) -> MatrixType { return solve(c, B); }, - nb::arg("B"), + "B"_a, "Returns the solution X of A X = B using the current " "decomposition of A where B is a right hand side matrix.") diff --git a/include/nanoeigenpy/decompositions/llt.hpp b/include/nanoeigenpy/decompositions/llt.hpp index c46f449..2a9d24c 100644 --- a/include/nanoeigenpy/decompositions/llt.hpp +++ b/include/nanoeigenpy/decompositions/llt.hpp @@ -8,6 +8,7 @@ namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template MatrixOrVector solve(const Eigen::LLT &c, @@ -16,7 +17,7 @@ MatrixOrVector solve(const Eigen::LLT &c, } template -void exposeLLTSolver(nb::module_ m, const char *name) { +void exposeLLT(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Chol = Eigen::LLT; using Scalar = typename MatrixType::Scalar; @@ -41,18 +42,18 @@ void exposeLLTSolver(nb::module_ m, const char *name) { "problems with hermitian matrices.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("size"), + .def(nb::init(), "size"_a, "Default constructor with memory preallocation.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LLT factorization from a given matrix.") .def(EigenBaseVisitor()) .def( - "matrixL", [](Chol const &c) -> MatrixType { return c.matrixL(); }, + "matrixL", [](const Chol &c) -> MatrixType { return c.matrixL(); }, "Returns the lower triangular matrix L.") .def( - "matrixU", [](Chol const &c) -> MatrixType { return c.matrixU(); }, + "matrixU", [](const Chol &c) -> MatrixType { return c.matrixU(); }, "Returns the upper triangular matrix U.") .def("matrixLLT", &Chol::matrixLLT, "Returns the LLT decomposition matrix made of the lower matrix " @@ -62,19 +63,18 @@ void exposeLLTSolver(nb::module_ m, const char *name) { #if EIGEN_VERSION_AT_LEAST(3, 3, 90) .def( "rankUpdate", - [](Chol &c, VectorType const &w, Scalar sigma) -> Chol & { + [](Chol &c, const VectorType &w, Scalar sigma) -> Chol & { return c.rankUpdate(w, sigma); }, - "If LL^* = A, then it becomes A + sigma * v v^*", nb::arg("w"), - nb::arg("sigma"), nb::rv_policy::reference) + "If LL^* = A, then it becomes A + sigma * v v^*", "w"_a, "sigma"_a, + nb::rv_policy::reference) #else .def( "rankUpdate", - [](Chol &c, VectorType const &w, Scalar sigma) -> Chol & { + [](Chol &c, const VectorType &w, Scalar sigma) -> Chol & { return c.rankUpdate(w, sigma); }, - "If LL^* = A, then it becomes A + sigma * v v^*", nb::arg("w"), - nb::arg("sigma")) + "If LL^* = A, then it becomes A + sigma * v v^*", "w"_a, "sigma"_a) #endif .def("adjoint", &Chol::adjoint, @@ -84,10 +84,10 @@ void exposeLLTSolver(nb::module_ m, const char *name) { .def( "compute", - [](Chol &c, MatrixType const &matrix) -> Chol & { + [](Chol &c, const MatrixType &matrix) -> Chol & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the LDLT of given matrix.", + "matrix"_a, "Computes the LDLT of given matrix.", nb::rv_policy::reference) .def("info", &Chol::info, "NumericalIssue if the input contains INF or NaN values or " @@ -104,18 +104,18 @@ void exposeLLTSolver(nb::module_ m, const char *name) { .def( "solve", - [](Chol const &c, VectorType const &b) -> VectorType { + [](const Chol &c, const VectorType &b) -> VectorType { return solve(c, b); }, - nb::arg("b"), + "b"_a, "Returns the solution x of A x = b using the current " "decomposition of A.") .def( "solve", - [](Chol const &c, MatrixType const &B) -> MatrixType { + [](const Chol &c, const MatrixType &B) -> MatrixType { return solve(c, B); }, - nb::arg("B"), + "B"_a, "Returns the solution X of A X = B using the current " "decomposition of A where B is a right hand side matrix.") diff --git a/include/nanoeigenpy/decompositions/partial-piv-lu.hpp b/include/nanoeigenpy/decompositions/partial-piv-lu.hpp new file mode 100644 index 0000000..6978875 --- /dev/null +++ b/include/nanoeigenpy/decompositions/partial-piv-lu.hpp @@ -0,0 +1,121 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +MatrixOrVector solve(const Eigen::PartialPivLU &c, + const MatrixOrVector &vec) { + return c.solve(vec); +} + +template +void exposePartialPivLU(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::PartialPivLU; + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + nb::class_( + m, name, + "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().") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "Constructs a LU factorization from a given matrix.") + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes the LU of given matrix.", + nb::rv_policy::reference) + + .def("matrixLU", &Solver::matrixLU, + "Returns the LU decomposition matrix: the upper-triangular part is " + "U, the " + "unit-lower-triangular part is L (at least for square matrices; in " + "the non-square " + "case, special care is needed, see the documentation of class " + "FullPivLU).", + nb::rv_policy::reference_internal) + + .def("permutationP", &Solver::permutationP, + "Returns the permutation matrix P in the decomposition A = P^{-1} L " + "U Q^{-1}.", + nb::rv_policy::reference_internal) + + .def("rcond", &Solver::rcond, + "Returns an estimate of the reciprocal condition number of the " + "matrix.") + .def( + "inverse", [](const Solver &c) -> MatrixType { return c.inverse(); }, + "Returns the inverse of the matrix associated with the LU " + "decomposition.") + .def("determinant", &Solver::determinant, + "Returns the determinant of the underlying matrix from the " + "current factorization.") + .def("reconstructedMatrix", &Solver::reconstructedMatrix, + "Returns the matrix represented by the decomposition," + "i.e., it returns the product: P^{-1} L U." + "This function is provided for debug purpose.") + + .def("rows", &Solver::rows, "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, "Returns the number of cols of the matrix.") + + .def( + "solve", + [](const Solver &c, const VectorType &b) -> VectorType { + return solve(c, b); + }, + "b"_a, + "Returns the solution x of A x = b using the current " + "decomposition of A.") + .def( + "solve", + [](const Solver &c, const MatrixType &B) -> MatrixType { + return solve(c, B); + }, + "B"_a, + "Returns the solution X of A X = B using the current " + "decomposition of A where B is a right hand side matrix.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/permutation-matrix.hpp b/include/nanoeigenpy/decompositions/permutation-matrix.hpp index 8cd0924..85d0514 100644 --- a/include/nanoeigenpy/decompositions/permutation-matrix.hpp +++ b/include/nanoeigenpy/decompositions/permutation-matrix.hpp @@ -8,6 +8,7 @@ namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template @@ -27,9 +28,9 @@ void exposePermutationMatrix(nb::module_ m, const char *name) { "This class represents a permutation matrix, " "internally stored as a vector of integers.") - .def(nb::init(), nb::arg("size"), + .def(nb::init(), "size"_a, "Default constructor with memory preallocation.") - .def(nb::init(), nb::arg("indices"), + .def(nb::init(), "indices"_a, "The indices array has the meaning that the permutations sends " "each integer i to indices[i].\n" "It is your responsibility to check that the indices array that " @@ -45,12 +46,10 @@ void exposePermutationMatrix(nb::module_ m, const char *name) { "The stored array representing the permutation.") .def("applyTranspositionOnTheLeft", - &PermutationMatrix::applyTranspositionOnTheLeft, nb::arg("i"), - nb::arg("j"), + &PermutationMatrix::applyTranspositionOnTheLeft, "i"_a, "j"_a, "Multiplies self by the transposition (ij) on the left.") .def("applyTranspositionOnTheRight", - &PermutationMatrix::applyTranspositionOnTheRight, nb::arg("i"), - nb::arg("j"), + &PermutationMatrix::applyTranspositionOnTheRight, "i"_a, "j"_a, "Multiplies self by the transposition (ij) on the right.") .def( @@ -61,7 +60,7 @@ void exposePermutationMatrix(nb::module_ m, const char *name) { [](PermutationMatrix &self, Eigen::DenseIndex size) { self.setIdentity(size); }, - nb::arg("size"), + "size"_a, "Sets self to be the identity permutation matrix of given size.") .def("toDenseMatrix", &PermutationMatrix::toDenseMatrix, @@ -70,18 +69,18 @@ void exposePermutationMatrix(nb::module_ m, const char *name) { .def( "transpose", - [](PermutationMatrix const &self) -> PermutationMatrix { + [](const PermutationMatrix &self) -> PermutationMatrix { return self.transpose(); }, "Returns the tranpose permutation matrix.") .def( "inverse", - [](PermutationMatrix const &self) -> PermutationMatrix { + [](const PermutationMatrix &self) -> PermutationMatrix { return self.inverse(); }, "Returns the inverse permutation matrix.") - .def("resize", &PermutationMatrix::resize, nb::arg("size"), + .def("resize", &PermutationMatrix::resize, "size"_a, "Resizes to given size.") .def(nb::self * nb::self) diff --git a/include/nanoeigenpy/decompositions/real-qz.hpp b/include/nanoeigenpy/decompositions/real-qz.hpp new file mode 100644 index 0000000..fb2b74d --- /dev/null +++ b/include/nanoeigenpy/decompositions/real-qz.hpp @@ -0,0 +1,69 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeRealQZ(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::RealQZ; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Real QZ decomposition") + + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "A"_a, + "B"_a, "computeQZ"_a = true, + "Constructor; computes real QZ decomposition of given matrices.") + + .def("matrixQ", &Solver::matrixQ, + "Returns matrix Q in the QZ decomposition.", + nb::rv_policy::reference_internal) + .def("matrixZ", &Solver::matrixZ, + "Returns matrix Z in the QZ decomposition.", + nb::rv_policy::reference_internal) + .def("matrixS", &Solver::matrixS, + "Returns matrix S in the QZ decomposition.", + nb::rv_policy::reference_internal) + .def("matrixT", &Solver::matrixT, + "Returns matrix T in the QZ decomposition.", + nb::rv_policy::reference_internal) + + .def( + "compute", + [](Solver &c, const MatrixType &A, const MatrixType &B) -> Solver & { + return c.compute(A, B); + }, + "A"_a, "B"_a, "Computes QZ decomposition of given matrix. ", + nb::rv_policy::reference) + + .def( + "compute", + [](Solver &c, const MatrixType &A, const MatrixType &B, + bool computeQZ) -> Solver & { return c.compute(A, B, computeQZ); }, + "A"_a, "B"_a, "computeQZ"_a, + "Computes QZ decomposition of given matrix. ", + nb::rv_policy::reference) + + .def("info", &Solver::info, + "Reports whether previous computation was successful.") + + .def("iterations", &Solver::iterations, + "Returns number of performed QR-like iterations.") + .def("setMaxIterations", &Solver::setMaxIterations, + "Sets the maximum number of iterations allowed.", + nb::rv_policy::reference) + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/real-schur.hpp b/include/nanoeigenpy/decompositions/real-schur.hpp new file mode 100644 index 0000000..d4062c2 --- /dev/null +++ b/include/nanoeigenpy/decompositions/real-schur.hpp @@ -0,0 +1,74 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeRealSchur(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::RealSchur; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Real Schur decomposition") + + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "computeU"_a = true, + "Constructor; computes real Schur decomposition of given matrices.") + + .def("matrixU", &Solver::matrixU, + "Returns the orthogonal matrix in the Schur decomposition. ", + nb::rv_policy::reference_internal) + .def("matrixT", &Solver::matrixT, + "Returns the quasi-triangular matrix in the Schur decomposition.", + nb::rv_policy::reference_internal) + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes Schur decomposition of given matrix.", + nb::rv_policy::reference) + + .def( + "compute", + [](Solver &c, const MatrixType &matrix, bool computeU) -> Solver & { + return c.compute(matrix, computeU); + }, + "matrix"_a, "computeU"_a, + "Computes Schur decomposition of given matrix.", + nb::rv_policy::reference) + + .def( + "computeFromHessenberg", + [](Solver &c, const MatrixType &matrixH, const MatrixType &matrixQ, + bool computeU) -> Solver & { + return c.computeFromHessenberg(matrixH, matrixQ, computeU); + }, + "matrixH"_a, "matrixQ"_a, "computeU"_a, + "Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T", + nb::rv_policy::reference) + + .def("info", &Solver::info, + "Reports whether previous computation was successful.") + + .def("setMaxIterations", &Solver::setMaxIterations, + "Sets the maximum number of iterations allowed.", + nb::rv_policy::reference) + .def("getMaxIterations", &Solver::getMaxIterations, + "Returns the maximum number of iterations.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/self-adjoint-eigen-solver.hpp b/include/nanoeigenpy/decompositions/self-adjoint-eigen-solver.hpp index 29cce61..3abff0b 100644 --- a/include/nanoeigenpy/decompositions/self-adjoint-eigen-solver.hpp +++ b/include/nanoeigenpy/decompositions/self-adjoint-eigen-solver.hpp @@ -3,17 +3,18 @@ #pragma once #include "nanoeigenpy/fwd.hpp" -#include "nanoeigenpy/eigen-base.hpp" -#include #include namespace nanoeigenpy { namespace nb = nanobind; +using namespace nb::literals; template void exposeSelfAdjointEigenSolver(nb::module_ m, const char *name) { using MatrixType = _MatrixType; using Solver = Eigen::SelfAdjointEigenSolver; + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Matrix; if (check_registration_alias(m)) { return; @@ -21,50 +22,56 @@ void exposeSelfAdjointEigenSolver(nb::module_ m, const char *name) { nb::class_(m, name, "Self adjoint Eigen Solver") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("size"), + .def(nb::init(), "size"_a, "Default constructor with memory preallocation.") .def(nb::init(), - nb::arg("matrix"), nb::arg("options") = Eigen::ComputeEigenvectors, + "matrix"_a, "options"_a = Eigen::ComputeEigenvectors, "Computes eigendecomposition of given matrix") - .def("eigenvalues", &Solver::eigenvalues, - "Returns the eigenvalues of given matrix.", - nb::rv_policy::reference_internal) - .def("eigenvectors", &Solver::eigenvectors, - "Returns the eigenvectors of given matrix.", - nb::rv_policy::reference_internal) + .def( + "eigenvalues", + [](const Solver &c) -> const VectorType & { return c.eigenvalues(); }, + "Returns the eigenvalues of given matrix.", + nb::rv_policy::reference_internal) + .def( + "eigenvectors", + [](const Solver &c) -> const MatrixType & { + return c.eigenvectors(); + }, + "Returns the eigenvectors of given matrix.", + nb::rv_policy::reference_internal) .def( "compute", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.compute(matrix); }, - nb::arg("matrix"), "Computes the eigendecomposition of given matrix.", + "matrix"_a, "Computes the eigendecomposition of given matrix.", nb::rv_policy::reference) .def( "compute", - [](Solver &c, MatrixType const &matrix, int options) -> Solver & { + [](Solver &c, const MatrixType &matrix, int options) -> Solver & { return c.compute(matrix, options); }, - nb::arg("matrix"), nb::arg("options"), + "matrix"_a, "options"_a, "Computes the eigendecomposition of given matrix.", nb::rv_policy::reference) .def( "computeDirect", - [](Solver &c, MatrixType const &matrix) -> Solver & { + [](Solver &c, const MatrixType &matrix) -> Solver & { return c.computeDirect(matrix); }, - nb::arg("matrix"), + "matrix"_a, "Computes eigendecomposition of given matrix using a closed-form " "algorithm.", nb::rv_policy::reference) .def( "computeDirect", - [](Solver &c, MatrixType const &matrix, int options) -> Solver & { + [](Solver &c, const MatrixType &matrix, int options) -> Solver & { return c.computeDirect(matrix, options); }, - nb::arg("matrix"), nb::arg("options"), + "matrix"_a, "options"_a, "Computes eigendecomposition of given matrix using a closed-form " "algorithm.", nb::rv_policy::reference) diff --git a/include/nanoeigenpy/decompositions/sparse/accelerate/accelerate.hpp b/include/nanoeigenpy/decompositions/sparse/accelerate/accelerate.hpp index 61d244c..d86bc0f 100644 --- a/include/nanoeigenpy/decompositions/sparse/accelerate/accelerate.hpp +++ b/include/nanoeigenpy/decompositions/sparse/accelerate/accelerate.hpp @@ -39,10 +39,10 @@ struct AccelerateImplVisitor .def( "compute", - [](Solver& c, MatrixType const& matrix) { + [](Solver& c, const MatrixType& matrix) { return c.compute(matrix); }, - nb::arg("matrix"), + "matrix"_a, "Computes the sparse Cholesky decomposition of a given matrix.", nb::rv_policy::reference) diff --git a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-base.hpp b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-base.hpp index bf61fd0..4878114 100644 --- a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-base.hpp +++ b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-base.hpp @@ -7,6 +7,7 @@ #include namespace nanoeigenpy { +using namespace nb::literals; struct CholmodBaseVisitor : nb::def_visitor { template @@ -26,10 +27,10 @@ struct CholmodBaseVisitor : nb::def_visitor { .def( "compute", - [](Solver &self, MatrixType const &matrix) -> decltype(auto) { + [](Solver &self, const MatrixType &matrix) -> decltype(auto) { return self.compute(matrix); }, - nb::arg("matrix"), + "matrix"_a, "Computes the sparse Cholesky decomposition of a given matrix.", nb::rv_policy::reference) @@ -37,7 +38,7 @@ struct CholmodBaseVisitor : nb::def_visitor { "Returns the determinant of the underlying matrix from the " "current factorization.") - .def("factorize", &Solver::factorize, nb::arg("matrix"), + .def("factorize", &Solver::factorize, "matrix"_a, "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" @@ -51,7 +52,7 @@ struct CholmodBaseVisitor : nb::def_visitor { "NumericalIssue if the input contains INF or NaN values or " "overflow occured. Returns Success otherwise.") - .def("setShift", &Solver::setShift, nb::arg("offset"), + .def("setShift", &Solver::setShift, "offset"_a, "Sets the shift parameters that will be used to adjust the " "diagonal coefficients during the numerical factorization.\n" "During the numerical factorization, the diagonal coefficients " diff --git a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-decomposition.hpp b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-decomposition.hpp index d9c8713..4032996 100644 --- a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-decomposition.hpp +++ b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-decomposition.hpp @@ -17,7 +17,7 @@ struct CholmodDecompositionVisitor "Eigen::SparseSolverBase"); cl.def(CholmodBaseVisitor()) - .def("setMode", &Solver::setMode, nb::arg("mode"), + .def("setMode", &Solver::setMode, "mode"_a, "Set the mode for the Cholesky decomposition."); } }; diff --git a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-ldlt.hpp b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-ldlt.hpp index f070b18..32e3a0b 100644 --- a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-ldlt.hpp +++ b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-ldlt.hpp @@ -30,7 +30,7 @@ void exposeCholmodSimplicialLDLT(nb::module_ m, const char *name) { "The vectors or matrices X and B can be either dense or sparse.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LDLT factorization from a given matrix.") .def(CholmodBaseVisitor()); diff --git a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-llt.hpp b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-llt.hpp index 7a40dac..ba56ecf 100644 --- a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-llt.hpp +++ b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-simplicial-llt.hpp @@ -31,7 +31,7 @@ void exposeCholmodSimplicialLLT(nb::module_ m, const char *name) { "The vectors or matrices X and B can be either dense or sparse.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LDLT factorization from a given matrix.") .def(CholmodBaseVisitor()); diff --git a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-supernodal-llt.hpp b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-supernodal-llt.hpp index a1a7666..385a4da 100644 --- a/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-supernodal-llt.hpp +++ b/include/nanoeigenpy/decompositions/sparse/cholmod/cholmod-supernodal-llt.hpp @@ -28,7 +28,7 @@ void exposeCholmodSupernodalLLT(nb::module_ m, const char *name) { "vectors or matrices X and B can be either dense or sparse.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LDLT factorization from a given matrix.") .def(CholmodBaseVisitor()); diff --git a/include/nanoeigenpy/decompositions/sparse/simplicial-cholesky.hpp b/include/nanoeigenpy/decompositions/sparse/simplicial-cholesky.hpp index b38266f..47f3c4d 100644 --- a/include/nanoeigenpy/decompositions/sparse/simplicial-cholesky.hpp +++ b/include/nanoeigenpy/decompositions/sparse/simplicial-cholesky.hpp @@ -7,6 +7,7 @@ #include namespace nanoeigenpy { +using namespace nb::literals; struct SimplicialCholeskyVisitor : nb::def_visitor { template @@ -26,19 +27,19 @@ struct SimplicialCholeskyVisitor : nb::def_visitor { .def( "matrixL", - [](Solver const &self) -> MatrixType { return self.matrixL(); }, + [](const Solver &self) -> MatrixType { return self.matrixL(); }, "Returns the lower triangular matrix L.") .def( "matrixU", - [](Solver const &self) -> MatrixType { return self.matrixU(); }, + [](const Solver &self) -> MatrixType { return self.matrixU(); }, "Returns the upper triangular matrix U.") .def( "compute", - [](Solver &self, MatrixType const &matrix) -> decltype(auto) { + [](Solver &self, const MatrixType &matrix) -> decltype(auto) { return self.compute(matrix); }, - nb::arg("matrix"), + "matrix"_a, "Computes the sparse Cholesky decomposition of a given matrix.", nb::rv_policy::reference) @@ -46,7 +47,7 @@ struct SimplicialCholeskyVisitor : nb::def_visitor { "Returns the determinant of the underlying matrix from the " "current factorization.") - .def("factorize", &Solver::factorize, nb::arg("matrix"), + .def("factorize", &Solver::factorize, "matrix"_a, "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" @@ -58,8 +59,8 @@ struct SimplicialCholeskyVisitor : nb::def_visitor { "NumericalIssue if the input contains INF or NaN values or " "overflow occured. Returns Success otherwise.") - .def("setShift", &Solver::setShift, nb::arg("offset"), - nb::arg("scale") = RealScalar(1), + .def("setShift", &Solver::setShift, "offset"_a, + "scale"_a = RealScalar(1), "Sets the shift parameters that will be used to adjust the " "diagonal coefficients during the numerical factorization.\n" "During the numerical factorization, the diagonal coefficients " diff --git a/include/nanoeigenpy/decompositions/sparse/simplicial-ldlt.hpp b/include/nanoeigenpy/decompositions/sparse/simplicial-ldlt.hpp index efb8342..4233826 100644 --- a/include/nanoeigenpy/decompositions/sparse/simplicial-ldlt.hpp +++ b/include/nanoeigenpy/decompositions/sparse/simplicial-ldlt.hpp @@ -34,12 +34,12 @@ void exposeSimplicialLDLT(nb::module_ m, const char *name) { "P^-1.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LDLT factorization from a given matrix.") .def( "vectorD", - [](Solver const &self) -> DenseVectorXs { return self.vectorD(); }, + [](const Solver &self) -> DenseVectorXs { return self.vectorD(); }, "Returns the diagonal vector D.") .def(SimplicialCholeskyVisitor()) diff --git a/include/nanoeigenpy/decompositions/sparse/simplicial-llt.hpp b/include/nanoeigenpy/decompositions/sparse/simplicial-llt.hpp index eb7e51a..d2093f7 100644 --- a/include/nanoeigenpy/decompositions/sparse/simplicial-llt.hpp +++ b/include/nanoeigenpy/decompositions/sparse/simplicial-llt.hpp @@ -30,7 +30,7 @@ void exposeSimplicialLLT(nb::module_ m, const char *name) { "P^-1.") .def(nb::init<>(), "Default constructor.") - .def(nb::init(), nb::arg("matrix"), + .def(nb::init(), "matrix"_a, "Constructs a LLT factorization from a given matrix.") .def(SimplicialCholeskyVisitor()) diff --git a/include/nanoeigenpy/decompositions/sparse/sparse-lu.hpp b/include/nanoeigenpy/decompositions/sparse/sparse-lu.hpp new file mode 100644 index 0000000..30d51b1 --- /dev/null +++ b/include/nanoeigenpy/decompositions/sparse/sparse-lu.hpp @@ -0,0 +1,181 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "nanoeigenpy/decompositions/sparse/sparse-solver-base.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +static void solveInPlace(const LTypeOrUType &self, + Eigen::Ref mat_vec) { + self.solveInPlace(mat_vec); +} + +template +void exposeMatrixL(nb::module_ m) { + using LType = Eigen::SparseLUMatrixLReturnType; + using Scalar = typename MappedSupernodalType::Scalar; + using VectorXs = Eigen::Matrix; + using MatrixXs = + Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + + nb::class_(m, "SparseLUMatrixLReturnType") + .def("rows", <ype::rows) + .def("cols", <ype::cols) + .def("solveInPlace", &solveInPlace, "X"_a) + .def("solveInPlace", &solveInPlace, "x"_a); +} + +template +void exposeMatrixU(nb::module_ m) { + using UType = Eigen::SparseLUMatrixUReturnType; + using Scalar = typename MatrixLType::Scalar; + using VectorXs = Eigen::Matrix; + using MatrixXs = + Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + + nb::class_(m, "SparseLUMatrixUReturnType") + .def("rows", &UType::rows) + .def("cols", &UType::cols) + .def("solveInPlace", &solveInPlace, "X"_a) + .def("solveInPlace", &solveInPlace, "x"_a); +} + +template > +void exposeSparseLU(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::SparseLU; + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + using StorageIndex = typename MatrixType::StorageIndex; + using SCMatrix = typename Solver::SCMatrix; + using MappedSparseMatrix = + typename Eigen::MappedSparseMatrix; + using LType = Eigen::SparseLUMatrixLReturnType; + using UType = Eigen::SparseLUMatrixUReturnType; + + if (check_registration_alias(m)) { + return; + } + + exposeMatrixL(m); + exposeMatrixU(m); + + nb::class_( + m, name, + "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.\n\n" + "An important parameter of this class is the ordering method. It is used " + "to reorder the " + "columns (and eventually the rows) of the matrix to reduce the number of " + "new elements that " + "are created during numerical factorization. The cheapest method " + "available is COLAMD. See " + "the OrderingMethods module for the list of built-in and external " + "ordering methods.") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "matrix"_a, + "Constructs a LU factorization from a given matrix.") + + .def(SparseSolverBaseVisitor()) + + .def("analyzePattern", &Solver::analyzePattern, + "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("factorize", &Solver::factorize, + "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("compute", &Solver::compute, + "Compute the symbolic and numeric factorization of the input sparse " + "matrix.\n\n" + "The input matrix should be in column-major storage.") + + .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("rows", &Solver::rows, "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, "Returns the number of cols of the matrix.") + + .def("isSymmetric", &Solver::isSymmetric, + "Indicate that the pattern of the input matrix is symmetric.") + + .def("rowsPermutation", &Solver::rowsPermutation, + "Returns a reference to the row matrix permutation " + "\f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$.", + nb::rv_policy::reference_internal) + .def("colsPermutation", &Solver::colsPermutation, + "Returns a reference to the column matrix permutation" + "\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$.", + nb::rv_policy::reference_internal) + + .def( + "setPivotThreshold", + [](Solver &self, const RealScalar &thresh) -> void { + return self.setPivotThreshold(thresh); + }, + "Set the threshold used for a diagonal entry to be an acceptable " + "pivot.") + + .def("info", &Solver::info, + "Reports whether previous computation was successful.") + .def("lastErrorMessage", &Solver::lastErrorMessage, + "A string describing the type of error") + + .def( + "absDeterminant", &Solver::absDeterminant, + "Returns the absolute value of the determinant of the matrix of which" + "*this is the QR decomposition.") + .def("logAbsDeterminant", &Solver::logAbsDeterminant, + "Returns the natural log of the absolute value of the determinant " + "of the matrix of which **this is the QR decomposition") + .def("signDeterminant", &Solver::signDeterminant, + "A number representing the sign of the determinant") + .def("determinant", &Solver::determinant, + "The determinant of the matrix.") + + .def("nnzL", &Solver::nnzL, "The number of non zero elements in L") + .def("nnzU", &Solver::nnzU, "The number of non zero elements in L") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/sparse/sparse-qr.hpp b/include/nanoeigenpy/decompositions/sparse/sparse-qr.hpp new file mode 100644 index 0000000..5758c1d --- /dev/null +++ b/include/nanoeigenpy/decompositions/sparse/sparse-qr.hpp @@ -0,0 +1,176 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "nanoeigenpy/decompositions/sparse/sparse-solver-base.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeMatrixQ(nb::module_ m) { + using Scalar = typename SparseQRType::Scalar; + using QType = Eigen::SparseQRMatrixQReturnType; + using QTransposeType = + Eigen::SparseQRMatrixQTransposeReturnType; + using VectorXd = Eigen::VectorXd; + using MatrixXd = Eigen::MatrixXd; + using QRMatrixType = typename SparseQRType::QRMatrixType; + + if (!check_registration_alias(m)) { + nb::class_(m, "SparseQRMatrixQTransposeReturnType") + .def(nb::init(), "qr"_a) + + .def( + "__matmul__", + [](QTransposeType& self, const MatrixXd& other) -> MatrixXd { + return MatrixXd(self * other); + }, + "other"_a) + + .def( + "__matmul__", + [](QTransposeType& self, const VectorXd& other) -> VectorXd { + return VectorXd(self * other); + }, + "other"_a); + } + + if (!check_registration_alias(m)) { + nb::class_(m, "SparseQRMatrixQReturnType") + .def(nb::init(), "qr"_a) + + .def("rows", &QType::rows) + .def("cols", &QType::cols) + + .def( + "__matmul__", + [](QType& self, const MatrixXd& other) -> MatrixXd { + return MatrixXd(self * other); + }, + "other"_a) + + .def( + "__matmul__", + [](QType& self, const VectorXd& other) -> VectorXd { + return VectorXd(self * other); + }, + "other"_a) + + .def("adjoint", + [](const QType& self) -> QTransposeType { return self.adjoint(); }) + + .def("transpose", [](const QType& self) -> QTransposeType { + return self.transpose(); + }); + } +} + +template > +void exposeSparseQR(nb::module_ m, const char* name) { + using MatrixType = _MatrixType; + using Ordering = _Ordering; + using Solver = Eigen::SparseQR; + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + using QRMatrixType = typename Solver::QRMatrixType; + using QType = Eigen::SparseQRMatrixQReturnType; + + if (check_registration_alias(m)) { + return; + } + + exposeMatrixQ(m); + + nb::class_( + m, name, + "Sparse left-looking QR factorization with numerical column pivoting. " + "\n\n" + "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. Use colsPermutation() to get it.\n\n" + "Q is the orthogonal matrix represented as products of Householder " + "reflectors. " + "Use matrixQ() to get an expression and matrixQ().adjoint() to get the " + "adjoint. Y" + "ou can then apply it to a vector.\n\n" + "R is the sparse triangular or trapezoidal matrix. The later occurs when " + "A is " + "rank-deficient. matrixR().topLeftCorner(rank(), rank()) always returns " + "a triangular " + "factor of full rank.") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "matrix"_a, + "Constructs a LU factorization from a given matrix.") + + .def(SparseSolverBaseVisitor()) + + .def("analyzePattern", &Solver::analyzePattern, + "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("factorize", &Solver::factorize, + "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("compute", &Solver::compute, + "Compute the symbolic and numeric factorization of the input sparse " + "matrix.\n\n" + "The input matrix should be in compressed mode " + "(see SparseMatrix::makeCompressed()).") + + .def( + "matrixQ", [](const Solver& self) -> QType { return self.matrixQ(); }, + "Returns an expression of the matrix Q as products of sparse " + "Householder reflectors.") + .def( + "matrixR", + [](const Solver& self) -> const QRMatrixType& { + return self.matrixR(); + }, + "Returns a const reference to the \b sparse upper triangular matrix " + "R of the QR factorization.", + nb::rv_policy::reference_internal) + + .def("rows", &Solver::rows, "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, "Returns the number of cols of the matrix.") + + .def("rank", &Solver::rank, + "Returns the number of non linearly dependent columns as determined " + "by the pivoting threshold.") + + .def("colsPermutation", &Solver::colsPermutation, + "Returns a reference to the column matrix permutation" + "\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$.", + nb::rv_policy::reference_internal) + + .def("info", &Solver::info, + "Reports whether previous computation was successful.") + .def("lastErrorMessage", &Solver::lastErrorMessage, + "A string describing the type of error") + + .def( + "setPivotThreshold", + [](Solver& self, const RealScalar& thresh) -> void { + return self.setPivotThreshold(thresh); + }, + "Set the threshold used for a diagonal entry to be an acceptable " + "pivot.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/sparse/sparse-solver-base.hpp b/include/nanoeigenpy/decompositions/sparse/sparse-solver-base.hpp index 549871e..36883fa 100644 --- a/include/nanoeigenpy/decompositions/sparse/sparse-solver-base.hpp +++ b/include/nanoeigenpy/decompositions/sparse/sparse-solver-base.hpp @@ -30,21 +30,21 @@ struct SparseSolverBaseVisitor : nb::def_visitor { cl.def( "solve", - [](Solver const &self, const Eigen::Ref &b) + [](const Solver &self, const Eigen::Ref &b) -> DenseVectorXs { return self.solve(b); }, "b"_a, "Returns the solution x of A x = b using the current decomposition " "of A, where b is a right hand side vector.") .def( "solve", - [](Solver const &self, const Eigen::Ref &B) + [](const Solver &self, const Eigen::Ref &B) -> DenseMatrixXs { return self.solve(B); }, "B"_a, "Returns the solution X of A X = B using the current decomposition " "of A where B is a right hand side matrix.") .def( "solve", - [](Solver const &self, const SparseMatrixType &B) + [](const Solver &self, const SparseMatrixType &B) -> SparseMatrixType { return self.solve(B); }, "B"_a, "Returns the solution X of A X = B using the current decomposition " diff --git a/include/nanoeigenpy/decompositions/svd-base.hpp b/include/nanoeigenpy/decompositions/svd-base.hpp new file mode 100644 index 0000000..680feee --- /dev/null +++ b/include/nanoeigenpy/decompositions/svd-base.hpp @@ -0,0 +1,83 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +struct SVDBaseVisitor : nb::def_visitor { + template + void execute(nb::class_ &cl) { + using MatrixType = typename Derived::MatrixType; + using RealScalar = typename MatrixType::RealScalar; + + using SVDBase = Eigen::SVDBase; + static_assert(std::is_base_of_v); + + cl.def(nb::init<>(), "Default constructor.") + + .def( + "matrixU", + [](const Derived &c) -> MatrixType { return c.matrixU(); }, + "Returns the U matrix") + .def( + "matrixV", + [](const Derived &c) -> MatrixType { return c.matrixV(); }, + "Returns the U matrix") + + .def("singularValues", &SVDBase::singularValues, + "For the SVD decomposition of a n-by-p matrix, letting " + "\a m be the minimum of \a n and \a p, the returned vector " + "has size \a m. Singular values are always sorted in decreasing " + "order.", + nb::rv_policy::reference_internal) + .def("nonzeroSingularValues", &SVDBase::nonzeroSingularValues, + "Returns the number of singular values that are not exactly 0.") + .def("rank", &SVDBase::rank, + "Returns the rank of the matrix of which *this is the SVD.") + + .def( + "setThreshold", + [](Derived &c, const RealScalar &threshold) { + return c.setThreshold(threshold); + }, + "threshold"_a, + "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 SVD decomposition " + "itself.\n\n" + "When it needs to get the threshold value, Eigen calls " + "threshold().", + nb::rv_policy::reference) + .def( + "setThreshold", + [](Derived &c) { return c.setThreshold(Eigen::Default); }, + "Allows to come back to the default behavior, letting Eigen use " + "its default formula for determining the threshold.", + nb::rv_policy::reference) + .def("threshold", &SVDBase::threshold, + "Returns the threshold that will be used by certain methods such " + "as rank().") + + .def("computeU", &SVDBase::computeU, + "Returns true if U (full or thin) is asked for in this SVD " + "decomposition") + .def("computeV", &SVDBase::computeV, + "Returns true if V (full or thin) is asked for in this SVD " + "decomposition") + + .def("rows", &SVDBase::rows, + "Returns the number of rows of the matrix.") + .def("cols", &SVDBase::cols, + "Returns the number of cols of the matrix.") + + .def("info", &SVDBase::info, + "Reports whether previous computation was successful."); + } +}; + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/decompositions/tridiagonalization.hpp b/include/nanoeigenpy/decompositions/tridiagonalization.hpp new file mode 100644 index 0000000..517bd0b --- /dev/null +++ b/include/nanoeigenpy/decompositions/tridiagonalization.hpp @@ -0,0 +1,60 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeTridiagonalization(nb::module_ m, const char *name) { + using MatrixType = _MatrixType; + using Solver = Eigen::Tridiagonalization; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Tridiagonalization") + + .def(nb::init(), "size"_a, + "Default constructor with memory preallocation.") + .def(nb::init(), "matrix"_a, + "Constructor; computes tridiagonal decomposition of given matrix.") + + .def( + "compute", + [](Solver &c, const MatrixType &matrix) -> Solver & { + return c.compute(matrix); + }, + "matrix"_a, "Computes tridiagonal decomposition of given matrix.", + nb::rv_policy::reference) + + .def("householderCoefficients", &Solver::householderCoefficients, + "Returns the Householder coefficients.") + + .def("packedMatrix", &Solver::packedMatrix, + "Returns the internal representation of the decomposition.", + nb::rv_policy::reference_internal) + + .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", &Solver::diagonal, + "Returns the diagonal of the tridiagonal matrix T in the " + "decomposition.") + .def("subDiagonal", &Solver::subDiagonal, + "Returns the subdiagonal of the tridiagonal matrix T in the " + "decomposition.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/fwd.hpp b/include/nanoeigenpy/fwd.hpp index d52821e..24b317c 100644 --- a/include/nanoeigenpy/fwd.hpp +++ b/include/nanoeigenpy/fwd.hpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include #if defined(__clang__) #define NANOEIGENPY_CLANG_COMPILER diff --git a/include/nanoeigenpy/geometry.hpp b/include/nanoeigenpy/geometry.hpp index 734d7ea..d9c7ee0 100644 --- a/include/nanoeigenpy/geometry.hpp +++ b/include/nanoeigenpy/geometry.hpp @@ -3,4 +3,10 @@ #pragma once #include "nanoeigenpy/geometry/angle-axis.hpp" +#include "nanoeigenpy/geometry/hyperplane.hpp" +#include "nanoeigenpy/geometry/parametrized-line.hpp" +#include "nanoeigenpy/geometry/rotation-2d.hpp" +#include "nanoeigenpy/geometry/scaling.hpp" +#include "nanoeigenpy/geometry/translation.hpp" #include "nanoeigenpy/geometry/quaternion.hpp" +#include "nanoeigenpy/geometry/jacobi-rotation.hpp" diff --git a/include/nanoeigenpy/geometry/angle-axis.hpp b/include/nanoeigenpy/geometry/angle-axis.hpp index 8278e80..4d3005a 100644 --- a/include/nanoeigenpy/geometry/angle-axis.hpp +++ b/include/nanoeigenpy/geometry/angle-axis.hpp @@ -88,7 +88,7 @@ void exposeAngleAxis(nb::module_ m, const char *name) { }, "other"_a, "Returns true if *this is approximately equal to other, " - "within the precision determined by prec.") + "within the default precision.") .def(nb::self * nb::self) .def(nb::self * Quaternion()) diff --git a/include/nanoeigenpy/geometry/hyperplane.hpp b/include/nanoeigenpy/geometry/hyperplane.hpp new file mode 100644 index 0000000..824a666 --- /dev/null +++ b/include/nanoeigenpy/geometry/hyperplane.hpp @@ -0,0 +1,192 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; + +template +bool isApprox( + const Eigen::Hyperplane &h, + const Eigen::Hyperplane &other, + const Scalar &prec = Eigen::NumTraits::dummy_precision()) { + return h.isApprox(other, prec); +} + +template +void exposeHyperplane(nb::module_ m, const char *name) { + using namespace nb::literals; + using Hyperplane = Eigen::Hyperplane; + using RealScalar = typename Eigen::NumTraits::Real; + using VectorType = Eigen::Matrix; + using Parameterized = Eigen::ParametrizedLine; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, + "A hyperplane is an affine subspace of " + "dimension n-1 in a space of dimension n.") + .def(nb::init<>(), "Default constructor") + .def(nb::init(), "copy"_a, "Copy constructor.") + .def(nb::init(), "dim"_a, + "Constructs a dynamic-size hyperplane with dim the dimension" + "of the ambient space.") + .def(nb::init(), "n"_a, "e"_a, + "Construct a plane from its normal \a n and a point \a e onto the " + "plane.") + .def(nb::init(), "n"_a, "d"_a, + "Constructs a plane from its normal n and distance to the origin d" + "such that the algebraic equation of the plane is \f$ n dot x + d = " + "0 \f$.") + .def(nb::init(), "parametrized"_a, + "Constructs a hyperplane passing through the parametrized line \a " + "parametrized." + "If the dimension of the ambient space is greater than 2, then " + "there isn't uniqueness," + "so an arbitrary choice is made.") + .def_static( + "Through", + [](const VectorType &p0, const VectorType &p1) -> Hyperplane { + return Hyperplane::Through(p0, p1); + }, + "p0"_a, "p1"_a, + "Constructs a hyperplane passing through the two points. " + "If the dimension of the ambient space is greater than 2, then " + "there isn't uniqueness, so an arbitrary choice is made.") + .def_static( + "Through", + [](const VectorType &p0, const VectorType &p1, + const VectorType &p2) -> Hyperplane { + if (p0.size() != 3 || p1.size() != 3 || p2.size() != 3) { + throw std::invalid_argument( + "Through with 3 points requires 3D vectors"); + } + + Hyperplane result(p0.size()); + VectorType v0 = p2 - p0; + VectorType v1 = p1 - p0; + + VectorType normal(3); + normal[0] = v0[1] * v1[2] - v0[2] * v1[1]; + normal[1] = v0[2] * v1[0] - v0[0] * v1[2]; + normal[2] = v0[0] * v1[1] - v0[1] * v1[0]; + + RealScalar norm = normal.norm(); + if (norm <= v0.norm() * v1.norm() * + Eigen::NumTraits::epsilon()) { + Eigen::Matrix m; + m.row(0) = v0.transpose(); + m.row(1) = v1.transpose(); + Eigen::JacobiSVD> svd( + m, Eigen::ComputeFullV); + result.normal() = svd.matrixV().col(2); + } else { + result.normal() = normal / norm; + } + + result.offset() = -p0.dot(result.normal()); + return result; + }, + "p0"_a, "p1"_a, "p2"_a, + "Constructs a hyperplane passing through the three points. " + "The dimension of the ambient space is required to be exactly 3.") + + .def("dim", &Hyperplane::dim, + "Returns the dimension in which the plane holds.") + .def("normalize", &Hyperplane::normalize, "Normalizes *this.") + + .def("signedDistance", &Hyperplane::signedDistance, "p"_a, + "Returns the signed distance between the plane *this and a point p.") + .def("absDistance", &Hyperplane::absDistance, "p"_a, + "Returns the absolute distance between the plane *this and a point " + "p.") + .def("projection", &Hyperplane::projection, "p"_a, + "Returns the projection of a point \a p onto the plane *this.") + + .def( + "normal", + [](const Hyperplane &self) -> VectorType { + return VectorType(self.normal()); + }, + "Returns a constant reference to the unit normal vector of the " + "plane, " + "which corresponds to the linear part of the implicit equation.") + .def( + "offset", + [](const Hyperplane &self) -> const Scalar & { + return self.offset(); + }, + "Returns the distance to the origin, which is also the constant " + "term of the implicit equation.", + nb::rv_policy::reference_internal) + .def( + "coeffs", [](const Hyperplane &self) { return self.coeffs(); }, + "Returns a constant reference to the coefficients c_i of the plane " + "equation: " + "\f$ c_0*x_0 + ... + c_{d-1}*x_{d-1} + c_d = 0 \f$.", + nb::rv_policy::reference_internal) + + .def( + "intersection", + [](const Hyperplane &self, const Hyperplane &other) -> VectorType { + if (self.dim() != 2 || other.dim() != 2) { + throw std::invalid_argument( + "intersection requires 2D hyperplanes"); + } + + Scalar det = self.coeffs().coeff(0) * other.coeffs().coeff(1) - + self.coeffs().coeff(1) * other.coeffs().coeff(0); + + if (Eigen::internal::isMuchSmallerThan(det, Scalar(1))) { + if (Eigen::numext::abs(self.coeffs().coeff(1)) > + Eigen::numext::abs(self.coeffs().coeff(0))) { + VectorType result(2); + result[0] = self.coeffs().coeff(1); + result[1] = -self.coeffs().coeff(2) / self.coeffs().coeff(1) - + self.coeffs().coeff(0); + return result; + } else { + VectorType result(2); + result[0] = -self.coeffs().coeff(2) / self.coeffs().coeff(0) - + self.coeffs().coeff(1); + result[1] = self.coeffs().coeff(0); + return result; + } + } else { + Scalar invdet = Scalar(1) / det; + VectorType result(2); + result[0] = + invdet * (self.coeffs().coeff(1) * other.coeffs().coeff(2) - + other.coeffs().coeff(1) * self.coeffs().coeff(2)); + result[1] = + invdet * (other.coeffs().coeff(0) * self.coeffs().coeff(2) - + self.coeffs().coeff(0) * other.coeffs().coeff(2)); + return result; + } + }, + "other"_a, "Returns the intersection of *this with \a other.") + + .def( + "isApprox", + [](const Hyperplane &aa, const Hyperplane &other, + const Scalar &prec) -> bool { return isApprox(aa, other, prec); }, + "other"_a, "prec"_a, + "Returns true if *this is approximately equal to other, " + "within the precision determined by prec.") + .def( + "isApprox", + [](const Hyperplane &aa, const Hyperplane &other) -> bool { + return isApprox(aa, other); + }, + "other"_a, + "Returns true if *this is approximately equal to other, " + "within the default precision.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/geometry/jacobi-rotation.hpp b/include/nanoeigenpy/geometry/jacobi-rotation.hpp new file mode 100644 index 0000000..bb5ae6d --- /dev/null +++ b/include/nanoeigenpy/geometry/jacobi-rotation.hpp @@ -0,0 +1,67 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeJacobiRotation(nb::module_ m, const char* name) { + using JacobiRotation = Eigen::JacobiRotation; + using RealScalar = typename Eigen::NumTraits::Real; + + if (check_registration_alias(m)) { + return; + } + nb::class_( + m, name, "This class represents a Jacobi or Givens rotation.") + .def(nb::init<>(), "Default constructor") + .def(nb::init(), "c"_a, "s"_a, + "Construct a planar rotation from a cosine-sine pair (c, s).") + + .def_prop_rw( + "c", [](const JacobiRotation& self) { return self.c(); }, + [](JacobiRotation& self, const Scalar& value) { self.c() = value; }) + .def_prop_rw( + "s", [](const JacobiRotation& self) { return self.s(); }, + [](JacobiRotation& self, const Scalar& value) { self.s() = value; }) + + .def("__mul__", &JacobiRotation::operator*, "other"_a, + "Concatenates two planar rotations") + + .def("transpose", &JacobiRotation::transpose) + .def("adjoint", &JacobiRotation::adjoint) + + .def( + "makeJacobi", + [](JacobiRotation& self, const RealScalar& x, const Scalar& y, + const RealScalar& z) { return self.makeJacobi(x, y, z); }, + "x"_a, "y"_a, "z"_a) + + .def( + "makeJacobi", + [](JacobiRotation& self, const Eigen::MatrixXd& m, Eigen::Index p, + Eigen::Index q) { return self.makeJacobi(m, p, q); }, + "matrix"_a, "p"_a, "q"_a) + + .def( + "makeGivens", + [](JacobiRotation& self, const Scalar& p, const Scalar& q) { + self.makeGivens(p, q, nullptr); + }, + "p"_a, "q"_a) + + .def( + "makeGivens", + [](JacobiRotation& self, const Scalar& p, const Scalar& q, + Scalar* r) { self.makeGivens(p, q, r); }, + "p"_a, "q"_a, "r"_a) + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/geometry/parametrized-line.hpp b/include/nanoeigenpy/geometry/parametrized-line.hpp new file mode 100644 index 0000000..3f9b5d3 --- /dev/null +++ b/include/nanoeigenpy/geometry/parametrized-line.hpp @@ -0,0 +1,139 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; + +template +bool isApprox( + const Eigen::ParametrizedLine &h, + const Eigen::ParametrizedLine &other, + const Scalar &prec = Eigen::NumTraits::dummy_precision()) { + return h.isApprox(other, prec); +} + +template +void exposeParametrizedLine(nb::module_ m, const char *name) { + using namespace nb::literals; + using ParametrizedLine = Eigen::ParametrizedLine; + using VectorType = Eigen::Matrix; + using Hyperplane = Eigen::Hyperplane; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Parametrized line.") + .def(nb::init<>(), "Default constructor") + .def(nb::init(), "copy"_a, "Copy constructor.") + .def(nb::init(), "dim"_a, + "Constructs a dynamic-size line with \a _dim the dimension" + "of the ambient space.") + .def(nb::init(), "origin"_a, + "direction"_a, + "Initializes a parametrized line of direction \a direction and " + "origin \a origin.") + .def(nb::init(), "copy"_a, "Copy constructor.") + .def( + "__init__", + [](ParametrizedLine *self, const Hyperplane &hyperplane) { + if (hyperplane.dim() != 2) { + throw std::invalid_argument( + "ParametrizedLine from Hyperplane requires ambient space " + "dimension 2"); + } + + VectorType normal = hyperplane.normal().eval(); + VectorType direction(2); + direction[0] = -normal[1]; + direction[1] = normal[0]; + + VectorType origin = -normal * hyperplane.offset(); + + new (self) ParametrizedLine(origin, direction); + }, + "hyperplane"_a, + "Constructs a parametrized line from a 2D hyperplane.") + .def_static( + "Through", + [](const VectorType &p0, const VectorType &p1) -> ParametrizedLine { + return ParametrizedLine::Through(p0, p1); + }, + "p0"_a, "p1"_a, + "Constructs a parametrized line going from \a p0 to \a p1.") + + .def("dim", &ParametrizedLine::dim, + "Returns the dimension in which the line holds.") + + .def( + "origin", + [](const ParametrizedLine &self) -> const VectorType & { + return self.origin(); + }, + nb::rv_policy::reference_internal, + "Returns the origin of the parametrized line.") + .def( + "direction", + [](const ParametrizedLine &self) -> const VectorType & { + return self.direction(); + }, + nb::rv_policy::reference_internal, + "Returns the direction of the parametrized line.") + + .def("squaredDistance", &ParametrizedLine::squaredDistance, "p"_a, + "Returns the squared distance of a point \a p to " + "its projection onto the line *this.") + .def("distance", &ParametrizedLine::distance, "p"_a, + "Returns the distance of a point \a p to " + "its projection onto the line *this.") + + .def("projection", &ParametrizedLine::projection, "p"_a, + "Returns the projection of a point \a p onto the line *this.") + .def("pointAt", &ParametrizedLine::pointAt, "t"_a) + + .def( + "intersectionParameter", + [](const ParametrizedLine &self, const Hyperplane &hyperplane) + -> Scalar { return self.intersectionParameter(hyperplane); }, + "hyperplane"_a, + "Returns the parameter value of the intersection between this line " + "and the given hyperplane.") + + .def( + "intersection", + [](const ParametrizedLine &self, const Hyperplane &hyperplane) + -> Scalar { return self.intersection(hyperplane); }, + "hyperplane"_a, + "Deprecated: use intersectionParameter(). Returns the parameter " + "value of the intersection.") + + .def( + "intersectionPoint", + [](const ParametrizedLine &self, const Hyperplane &hyperplane) + -> VectorType { return self.intersectionPoint(hyperplane); }, + "hyperplane"_a, + "Returns the point of the intersection between this line and the " + "given hyperplane.") + + .def( + "isApprox", + [](const ParametrizedLine &aa, const ParametrizedLine &other, + const Scalar &prec) -> bool { return isApprox(aa, other, prec); }, + "other"_a, "prec"_a, + "Returns true if *this is approximately equal to other, " + "within the precision determined by prec.") + .def( + "isApprox", + [](const ParametrizedLine &aa, const ParametrizedLine &other) + -> bool { return isApprox(aa, other); }, + "other"_a, + "Returns true if *this is approximately equal to other, " + "within the default precision.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/geometry/quaternion.hpp b/include/nanoeigenpy/geometry/quaternion.hpp index 17ae421..f5dd209 100644 --- a/include/nanoeigenpy/geometry/quaternion.hpp +++ b/include/nanoeigenpy/geometry/quaternion.hpp @@ -62,7 +62,6 @@ struct QuaternionVisitor : nb::def_visitor> { }, "u"_a, "v"_a, "Initialize from two vectors u and v.") - // Add property with x, y, z, w .def_prop_rw("x", &QuaternionVisitor::getCoeff<0>, &QuaternionVisitor::setCoeff<0>, "The x coefficient.") .def_prop_rw("y", &QuaternionVisitor::getCoeff<1>, @@ -162,7 +161,7 @@ struct QuaternionVisitor : nb::def_visitor> { [](Quaternion* self, const AngleAxisType& aa) -> Quaternion& { return (*self = aa); }, - nb::arg("aa"), nb::rv_policy::reference, + "aa"_a, nb::rv_policy::reference, "Set *this from an angle-axis and return a reference to self.") .def("__str__", &print) .def("__repr__", &print) diff --git a/include/nanoeigenpy/geometry/rotation-2d.hpp b/include/nanoeigenpy/geometry/rotation-2d.hpp new file mode 100644 index 0000000..207a110 --- /dev/null +++ b/include/nanoeigenpy/geometry/rotation-2d.hpp @@ -0,0 +1,121 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "detail/rotation-base.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; + +template +bool isApprox( + const Eigen::Rotation2D &r, const Eigen::Rotation2D &other, + const Scalar &prec = Eigen::NumTraits::dummy_precision()) { + return r.isApprox(other, prec); +} + +template +void exposeRotation2D(nb::module_ m, const char *name) { + using namespace nb::literals; + using Rotation2D = Eigen::Rotation2D; + using Vector2 = typename Rotation2D::Vector2; + using Matrix2 = typename Rotation2D::Matrix2; + + if (check_registration_alias(m)) { + return; + } + nb::class_( + m, name, "Represents a rotation/orientation in a 2 dimensional space.") + .def(nb::init<>(), "Default constructor") + .def(nb::init(), "a"_a, + "Construct a 2D counter clock wise rotation from the angle \a a in " + "radian.") + .def(nb::init(), "mat"_a, + "Construct a 2D rotation from a 2x2 rotation matrix \a mat.") + .def(nb::init(), "copy"_a, "Copy constructor.") + + .def_prop_rw( + "angle", [](const Rotation2D &r) -> Scalar { return r.angle(); }, + [](Rotation2D &r, Scalar a) { r.angle() = a; }, "The rotation angle.") + + .def("smallestPositiveAngle", &Rotation2D::smallestPositiveAngle, + "Returns the rotation angle in [0,2pi]") + .def("smallestAngle", &Rotation2D::smallestAngle, + "Returns the rotation angle in [-pi,pi]") + + .def(RotationBaseVisitor()) + + .def( + "fromRotationMatrix", + [](Rotation2D &r, const Matrix2 &mat) -> Rotation2D & { + return r.fromRotationMatrix(mat); + }, + "mat"_a, "Sets *this from a 2x2 rotation matrix", + nb::rv_policy::reference_internal) + .def("toRotationMatrix", &Rotation2D::toRotationMatrix) + + .def( + "slerp", + [](const Rotation2D &self, const Scalar t, const Rotation2D &other) + -> Rotation2D { return self.slerp(t, other); }, + "t"_a, "other"_a, + "Returns the spherical interpolation between *this and \a other using" + "parameter \a t. It is in fact equivalent to a linear interpolation.") + + .def_static("Identity", &Rotation2D::Identity) + + .def( + "isApprox", + [](const Rotation2D &r, const Rotation2D &other, + const Scalar &prec) -> bool { return isApprox(r, other, prec); }, + "other"_a, "prec"_a, + "Returns true if *this is approximately equal to other, " + "within the precision determined by prec.") + .def( + "isApprox", + [](const Rotation2D &r, const Rotation2D &other) -> bool { + return isApprox(r, other); + }, + "other"_a, + "Returns true if *this is approximately equal to other, " + "within the default precision.") + + .def( + "__mul__", + [](const Rotation2D &self, const Rotation2D &other) -> Rotation2D { + return self * other; + }, + "other"_a, "Concatenates two rotations") + .def( + "__imul__", + [](Rotation2D &self, const Rotation2D &other) -> Rotation2D & { + return self *= other; + }, + "other"_a, "Concatenates two rotations in-place", + nb::rv_policy::reference_internal) + .def( + "__mul__", + [](const Rotation2D &self, const Vector2 &vec) -> Vector2 { + return self * vec; + }, + "vec"_a, "Applies the rotation to a 2D vector") + + .def( + "__eq__", + [](const Rotation2D &self, const Rotation2D &other) -> bool { + return std::abs(self.angle() - other.angle()) < 1e-12; + }, + "other"_a, "Tests equality") + .def( + "__ne__", + [](const Rotation2D &self, const Rotation2D &other) -> bool { + return std::abs(self.angle() - other.angle()) >= 1e-12; + }, + "other"_a, "Tests inequality") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/geometry/scaling.hpp b/include/nanoeigenpy/geometry/scaling.hpp new file mode 100644 index 0000000..c92ccb7 --- /dev/null +++ b/include/nanoeigenpy/geometry/scaling.hpp @@ -0,0 +1,92 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "detail/rotation-base.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; + +template +bool isApprox( + const Eigen::UniformScaling& r, + const Eigen::UniformScaling& other, + const Scalar& prec = Eigen::NumTraits::dummy_precision()) { + return r.isApprox(other, prec); +} + +template +void exposeUniformScaling(nb::module_ m, const char* name) { + using namespace nb::literals; + using UniformScaling = Eigen::UniformScaling; + + if (check_registration_alias(m)) { + return; + } + nb::class_( + m, name, "Represents a generic uniform scaling transformation.") + .def(nb::init<>(), "Default constructor") + .def(nb::init(), "s"_a, + "Constructs and initialize a uniform scaling transformation") + .def(nb::init(), "copy"_a, "Copy constructor.") + + .def( + "factor", + [](const UniformScaling& self) -> const Scalar& { + return self.factor(); + }, + nb::rv_policy::reference_internal, "Returns the scaling factor") + .def("inverse", &UniformScaling::inverse, "Returns the inverse scaling") + + .def( + "isApprox", + [](const UniformScaling& r, const UniformScaling& other, + const Scalar& prec) -> bool { return isApprox(r, other, prec); }, + "other"_a, "prec"_a, + "Returns true if *this is approximately equal to other, " + "within the precision determined by prec.") + .def( + "isApprox", + [](const UniformScaling& r, const UniformScaling& other) -> bool { + return isApprox(r, other); + }, + "other"_a, + "Returns true if *this is approximately equal to other, " + "within the default precision.") + + .def( + "__mul__", + [](const UniformScaling& self, const UniformScaling& other) + -> UniformScaling { return self * other; }, + "other"_a, "Concatenates two uniform scalings") + + .def( + "__mul__", + [](const UniformScaling& self, const Eigen::MatrixXd& matrix) + -> Eigen::MatrixXd { return self * matrix; }, + "matrix"_a, "Multiplies uniform scaling with a matrix") + + .def( + "__mul__", + [](const UniformScaling& self, const Eigen::AngleAxis& r) + -> Eigen::Matrix { return self * r; }, + "rotation"_a, "Multiplies uniform scaling with AngleAxis rotation") + + .def( + "__mul__", + [](const UniformScaling& self, const Eigen::Quaternion& q) + -> Eigen::Matrix { return self * q; }, + "quaternion"_a, "Multiplies uniform scaling with quaternion") + + .def( + "__mul__", + [](const UniformScaling& self, const Eigen::Rotation2D& r) + -> Eigen::Matrix { return self * r; }, + "rotation2d"_a, "Multiplies uniform scaling with 2D rotation") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/geometry/translation.hpp b/include/nanoeigenpy/geometry/translation.hpp new file mode 100644 index 0000000..ec0f245 --- /dev/null +++ b/include/nanoeigenpy/geometry/translation.hpp @@ -0,0 +1,134 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include + +namespace nanoeigenpy { +namespace nb = nanobind; + +template +bool isApprox( + const Eigen::Translation& r, + const Eigen::Translation& other, + const Scalar& prec = Eigen::NumTraits::dummy_precision()) { + return r.isApprox(other, prec); +} + +template +void exposeTranslation(nb::module_ m, const char* name) { + using namespace nb::literals; + using VectorType = Eigen::Matrix; + using LinearMatrixType = + Eigen::Matrix; + using Translation = Eigen::Translation; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, "Represents a translation transformation.") + .def(nb::init<>(), "Default constructor") + .def(nb::init(), "sx"_a, "sy"_a) + .def(nb::init(), "sx"_a, + "sy"_a, "sz"_a) + .def(nb::init(), "vector"_a) + .def(nb::init(), "copy"_a, "Copy constructor.") + + .def_prop_rw( + "x", + [](const Translation& self) -> Scalar { + if (self.vector().size() < 1) { + throw std::out_of_range( + "Translation must have at least 1 dimension for x"); + } + return self.vector()[0]; + }, + [](Translation& self, const Scalar& value) { + if (self.vector().size() < 1) { + throw std::out_of_range( + "Translation must have at least 1 dimension for x"); + } + self.vector()[0] = value; + }, + "The x-translation") + .def_prop_rw( + "y", + [](const Translation& self) -> Scalar { + if (self.vector().size() < 2) { + throw std::out_of_range( + "Translation must have at least 2 dimensions for y"); + } + return self.vector()[1]; + }, + [](Translation& self, const Scalar& value) { + if (self.vector().size() < 2) { + throw std::out_of_range( + "Translation must have at least 2 dimensions for y"); + } + self.vector()[1] = value; + }, + "The y-translation") + .def_prop_rw( + "z", + [](const Translation& self) -> Scalar { + if (self.vector().size() < 3) { + throw std::out_of_range( + "Translation must have at least 3 dimensions for z"); + } + return self.vector()[2]; + }, + [](Translation& self, const Scalar& value) { + if (self.vector().size() < 3) { + throw std::out_of_range( + "Translation must have at least 3 dimensions for z"); + } + self.vector()[2] = value; + }, + "The z-translation") + + .def( + "vector", + [](const Translation& self) -> const VectorType& { + return self.vector(); + }, + nb::rv_policy::reference_internal, "Returns the translation vector") + + .def( + "translation", + [](const Translation& self) -> const VectorType& { + return self.translation(); + }, + nb::rv_policy::reference_internal, + "Returns the translation vector (alias for vector)") + + .def("inverse", &Translation::inverse, + "Returns the inverse translation (opposite)") + + .def( + "isApprox", + [](const Translation& r, const Translation& other, + const Scalar& prec) -> bool { return isApprox(r, other, prec); }, + "other"_a, "prec"_a, + "Returns true if *this is approximately equal to other, " + "within the precision determined by prec.") + .def( + "isApprox", + [](const Translation& r, const Translation& other) -> bool { + return isApprox(r, other); + }, + "other"_a, + "Returns true if *this is approximately equal to other, " + "within the default precision.") + + .def( + "__mul__", + [](const Translation& self, const Translation& other) -> Translation { + return self * other; + }, + "other"_a, "Concatenates two translations") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/id.hpp b/include/nanoeigenpy/id.hpp index d9355f8..027be19 100644 --- a/include/nanoeigenpy/id.hpp +++ b/include/nanoeigenpy/id.hpp @@ -12,7 +12,7 @@ struct IdVisitor : nb::def_visitor { void execute(nb::class_ &cl) { cl.def( "id", - [](C const &self) -> int64_t { + [](const C &self) -> int64_t { return reinterpret_cast(&self); }, "Returns the unique identity of an object.\n" diff --git a/include/nanoeigenpy/solvers.hpp b/include/nanoeigenpy/solvers.hpp index 8804c04..e1c5e05 100644 --- a/include/nanoeigenpy/solvers.hpp +++ b/include/nanoeigenpy/solvers.hpp @@ -10,3 +10,6 @@ #include "nanoeigenpy/solvers/least-squares-conjugate-gradient.hpp" #endif #include "nanoeigenpy/solvers/conjugate-gradient.hpp" +#include "nanoeigenpy/solvers/bicgstab.hpp" +#include "nanoeigenpy/solvers/incomplete-lut.hpp" +#include "nanoeigenpy/solvers/incomplete-cholesky.hpp" diff --git a/include/nanoeigenpy/solvers/basic-preconditioners.hpp b/include/nanoeigenpy/solvers/basic-preconditioners.hpp index 295eea3..45a9d8c 100644 --- a/include/nanoeigenpy/solvers/basic-preconditioners.hpp +++ b/include/nanoeigenpy/solvers/basic-preconditioners.hpp @@ -21,16 +21,15 @@ struct PreconditionerBaseVisitor .def(nb::init(), "A"_a) .def("info", &Preconditioner::info, "Returns success if the Preconditioner has been well initialized.") - .def("solve", &solve, nb::arg("b"), + .def("solve", &solve, "b"_a, "Returns the solution A * z = b where the preconditioner is an " "estimate of A^-1.") - .def("compute", &Preconditioner::template compute, - nb::arg("mat"), + .def("compute", &Preconditioner::template compute, "mat"_a, "Initialize the preconditioner from the matrix value.", nb::rv_policy::reference) .def("factorize", &Preconditioner::template factorize, - nb::arg("mat"), + "mat"_a, "Initialize the preconditioner from the matrix value, i.e " "factorize the mat given as input to approximate its inverse.", nb::rv_policy::reference); diff --git a/include/nanoeigenpy/solvers/bicgstab.hpp b/include/nanoeigenpy/solvers/bicgstab.hpp new file mode 100644 index 0000000..1959de8 --- /dev/null +++ b/include/nanoeigenpy/solvers/bicgstab.hpp @@ -0,0 +1,44 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" +#include "nanoeigenpy/solvers/iterative-solver-base.hpp" +#include +#include + +namespace nanoeigenpy { +namespace nb = nanobind; + +template +struct BiCGSTABVisitor : nb::def_visitor> { + using MatrixType = typename BiCGSTAB::MatrixType; + using CtorArg = nb::DMap; + + template + void execute(nb::class_& cl) { + using namespace nb::literals; + cl.def(nb::init<>(), "Default constructor.") + .def(nb::init(), "A"_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().") + .def(IterativeSolverVisitor()); + } + + static void expose(nb::module_& m, const char* name) { + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name) + .def(BiCGSTABVisitor()) + .def(IdVisitor()); + } +}; + +template +void exposeBiCGSTAB(nb::module_& m, const char* name) { + BiCGSTABVisitor::expose(m, name); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/solvers/incomplete-cholesky.hpp b/include/nanoeigenpy/solvers/incomplete-cholesky.hpp new file mode 100644 index 0000000..4952e4e --- /dev/null +++ b/include/nanoeigenpy/solvers/incomplete-cholesky.hpp @@ -0,0 +1,108 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeIncompleteCholesky(nb::module_ m, const char* name) { + using MatrixType = _MatrixType; + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename Eigen::NumTraits::Real; + using Solver = Eigen::IncompleteCholesky; + using Factortype = Eigen::SparseMatrix; + using VectorRx = Eigen::Matrix; + using PermutationType = + Eigen::PermutationMatrix; + + static constexpr int Options = + MatrixType::Options; // Options = Eigen::ColMajor + using DenseVectorXs = Eigen::Matrix; + using DenseMatrixXs = + Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + nb::class_(m, name, + "Modified Incomplete Cholesky with dual threshold.") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "matrix"_a, + "Constructs an incomplete cholesky factorization from a given " + "matrix.") + + .def("rows", &Solver::rows, "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, "Returns the number of cols of the matrix.") + + .def("info", &Solver::info, + "Reports whether previous computation was successful.") + + .def("setInitialShift", &Solver::setInitialShift, "shift"_a, + "Set the initial shift parameter.") + + .def( + "analyzePattern", + [](Solver& self, const MatrixType& amat) { + self.analyzePattern(amat); + }, + "matrix"_a) + .def( + "factorize", + [](Solver& self, const MatrixType& amat) { self.factorize(amat); }, + "matrix"_a) + .def( + "compute", + [](Solver& self, const MatrixType& amat) { self.compute(amat); }, + "matrix"_a) + + .def( + "matrixL", + [](const Solver& self) -> const Factortype& { + return self.matrixL(); + }, + nb::rv_policy::reference_internal) + .def( + "scalingS", + [](const Solver& self) -> const VectorRx& { return self.scalingS(); }, + nb::rv_policy::reference_internal) + .def( + "permutationP", + [](const Solver& self) -> const PermutationType& { + return self.permutationP(); + }, + nb::rv_policy::reference_internal) + + .def( + "solve", + [](const Solver& self, const Eigen::Ref& b) + -> DenseVectorXs { return self.solve(b); }, + "b"_a, + "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); }, + "B"_a, + "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()); + }, + "B"_a, + "Returns the solution X of A X = B using the current decomposition " + "of A where B is a right hand side matrix.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/solvers/incomplete-lut.hpp b/include/nanoeigenpy/solvers/incomplete-lut.hpp new file mode 100644 index 0000000..6cbc217 --- /dev/null +++ b/include/nanoeigenpy/solvers/incomplete-lut.hpp @@ -0,0 +1,113 @@ +/// Copyright 2025 INRIA + +#pragma once + +#include "nanoeigenpy/fwd.hpp" + +namespace nanoeigenpy { +namespace nb = nanobind; +using namespace nb::literals; + +template +void exposeIncompleteLUT(nb::module_ m, const char* name) { + using MatrixType = _MatrixType; + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename Eigen::NumTraits::Real; + using Solver = Eigen::IncompleteLUT; + + static constexpr int Options = + MatrixType::Options; // Options = Eigen::ColMajor + using DenseVectorXs = Eigen::Matrix; + using DenseMatrixXs = + Eigen::Matrix; + + if (check_registration_alias(m)) { + return; + } + nb::class_( + m, name, + "Incomplete LU factorization with dual-threshold strategy. " + "\n\n" + "This class follows the sparse solver concept .\n\n" + "During the numerical factorization, two dropping rules are used : " + "1) any element whose magnitude is less than some tolerance is dropped. " + "This tolerance is obtained by multiplying the input tolerance droptol " + "by the average magnitude of all the original elements in the current " + "row. 2) After the elimination of the row, only the fill largest " + "elements " + "in the L part and the fill largest elements in the U part are kept (in " + "addition to the diagonal element ). Note that fill is computed from the " + "input parameter fillfactor which is used the ratio to control the " + "fill_in " + "relatively to the initial number of nonzero elements. The two extreme " + "cases are when droptol=0 (to keep all the fill*2 largest elements) and " + "when fill=n/2 with droptol being different to zero.") + + .def(nb::init<>(), "Default constructor.") + .def(nb::init(), "matrix"_a, + "Constructs an incomplete LU factorization from a given matrix.") + .def( + "__init__", + [](Solver* self, const MatrixType& mat, + RealScalar droptol = Eigen::NumTraits::dummy_precision(), + int fillfactor = 10) { + new (self) Solver(mat, droptol, fillfactor); + }, + "matrix"_a, "droptol"_a = Eigen::NumTraits::dummy_precision(), + "fillfactor"_a = 10) + + .def("rows", &Solver::rows, "Returns the number of rows of the matrix.") + .def("cols", &Solver::cols, "Returns the number of cols of the matrix.") + + .def("info", &Solver::info, + "Reports whether previous computation was successful.") + + .def( + "analyzePattern", + [](Solver& self, const MatrixType& amat) { + self.analyzePattern(amat); + }, + "matrix"_a) + .def( + "factorize", + [](Solver& self, const MatrixType& amat) { self.factorize(amat); }, + "matrix"_a) + .def( + "compute", + [](Solver& self, const MatrixType& amat) -> Solver& { + return self.compute(amat); + }, + "matrix"_a, nb::rv_policy::reference) + + .def("setDroptol", &Solver::setDroptol) + .def("setFillfactor", &Solver::setFillfactor) + + .def( + "solve", + [](const Solver& self, const Eigen::Ref& b) + -> DenseVectorXs { return self.solve(b); }, + "b"_a, + "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); }, + "B"_a, + "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()); + }, + "B"_a, + "Returns the solution X of A X = B using the current decomposition " + "of A where B is a right hand side matrix.") + + .def(IdVisitor()); +} + +} // namespace nanoeigenpy diff --git a/include/nanoeigenpy/solvers/least-squares-conjugate-gradient.hpp b/include/nanoeigenpy/solvers/least-squares-conjugate-gradient.hpp index fdbf84f..fec82d4 100644 --- a/include/nanoeigenpy/solvers/least-squares-conjugate-gradient.hpp +++ b/include/nanoeigenpy/solvers/least-squares-conjugate-gradient.hpp @@ -24,7 +24,8 @@ struct LeastSquaresConjugateGradientVisitor "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()."); + "followed by a call to compute().") + .def(IterativeSolverVisitor()); } static void expose(nb::module_& m, const char* name) { @@ -32,7 +33,6 @@ struct LeastSquaresConjugateGradientVisitor return; } nb::class_(m, name) - .def(IterativeSolverVisitor()) .def(LeastSquaresConjugateGradientVisitor< LeastSquaresConjugateGradient>()) .def(IdVisitor()); diff --git a/include/nanoeigenpy/solvers/minres.hpp b/include/nanoeigenpy/solvers/minres.hpp index 6186b4d..bc32863 100644 --- a/include/nanoeigenpy/solvers/minres.hpp +++ b/include/nanoeigenpy/solvers/minres.hpp @@ -10,57 +10,33 @@ namespace nanoeigenpy { namespace nb = nanobind; -template -struct MINRESSolverVisitor : nb::def_visitor> { - using MatrixType = _MatrixType; - using Scalar = typename MatrixType::Scalar; - using RealScalar = typename MatrixType::RealScalar; - using VectorXs = - Eigen::Matrix; - using MatrixXs = Eigen::Matrix; - using Solver = Eigen::MINRES; - using CtorArg = nb::DMap; +template +struct MINRESVisitor : nb::def_visitor> { + using MatrixType = typename MINRES::MatrixType; + using CtorArg = nb::DMap; - public: template - void execute(nb::class_& cl) { + void execute(nb::class_& cl) { using namespace nb::literals; - cl.def(nb::init<>(), "Default constructor.") .def(nb::init(), "A"_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().") - .def(IterativeSolverVisitor()); + .def(IterativeSolverVisitor()); } static void expose(nb::module_& m, const char* name) { - if (check_registration_alias(m)) { + if (check_registration_alias(m)) { return; } - nb::class_( - m, name, - "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") - .def(MINRESSolverVisitor()) - .def(IdVisitor()); + nb::class_(m, name).def(MINRESVisitor()).def(IdVisitor()); } }; -template -void exposeMINRESSolver(nb::module_& m, const char* name) { - MINRESSolverVisitor<_MatrixType>::expose(m, name); +template +void exposeMINRES(nb::module_& m, const char* name) { + MINRESVisitor::expose(m, name); } } // namespace nanoeigenpy diff --git a/include/nanoeigenpy/utils/is-approx.hpp b/include/nanoeigenpy/utils/is-approx.hpp index a236cdd..ff4356b 100644 --- a/include/nanoeigenpy/utils/is-approx.hpp +++ b/include/nanoeigenpy/utils/is-approx.hpp @@ -9,6 +9,7 @@ #include namespace nanoeigenpy { +using namespace nb::literals; template EIGEN_DONT_INLINE bool is_approx( @@ -65,7 +66,7 @@ void exposeIsApprox(nb::module_ m) { [](const MatrixXs& mat1, const MatrixXs& mat2, RealScalar precision) { return is_approx(mat1, mat2, precision); }, - nb::arg("mat1"), nb::arg("mat2"), nb::arg("precision") = dummy_precision, + "mat1"_a, "mat2"_a, "precision"_a = dummy_precision, "Check if two dense matrices are approximately equal."); // is_approx for dense vectors @@ -74,7 +75,7 @@ void exposeIsApprox(nb::module_ m) { [](const VectorXs& vec1, const VectorXs& vec2, RealScalar precision) { return is_approx(vec1, vec2, precision); }, - nb::arg("vec1"), nb::arg("vec2"), nb::arg("precision") = dummy_precision, + "vec1"_a, "vec2"_a, "precision"_a = dummy_precision, "Check if two dense vectors are approximately equal."); } diff --git a/src/module.cpp b/src/module.cpp index cc478e1..3dfd59a 100644 --- a/src/module.cpp +++ b/src/module.cpp @@ -4,24 +4,67 @@ #include "nanoeigenpy/decompositions.hpp" #include "nanoeigenpy/geometry.hpp" -#include "nanoeigenpy/utils/is-approx.hpp" +#include "nanoeigenpy/solvers.hpp" #include "nanoeigenpy/constants.hpp" +#include "nanoeigenpy/utils/is-approx.hpp" #include "./internal.h" using namespace nanoeigenpy; -using Quaternion = Eigen::Quaternion; +// Matrix types using SparseMatrix = Eigen::SparseMatrix; +using Eigen::ColPivHouseholderQRPreconditioner; +using Eigen::FullPivHouseholderQRPreconditioner; +using Eigen::HouseholderQRPreconditioner; +using Eigen::JacobiSVD; +using Eigen::NoQRPreconditioner; + +using ColPivHhJacobiSVD = JacobiSVD; +using FullPivHhJacobiSVD = + JacobiSVD; +using HhJacobiSVD = JacobiSVD; +using NoPrecondJacobiSVD = JacobiSVD; + +using SparseQR = Eigen::SparseQR>; +using SparseLU = Eigen::SparseLU; +using SCMatrix = typename SparseLU::SCMatrix; +using StorageIndex = typename Matrix::StorageIndex; +using MappedSparseMatrix = + typename Eigen::MappedSparseMatrix; + +NB_MAKE_OPAQUE(ColPivHhJacobiSVD) +NB_MAKE_OPAQUE(FullPivHhJacobiSVD) +NB_MAKE_OPAQUE(HhJacobiSVD) +NB_MAKE_OPAQUE(NoPrecondJacobiSVD) + +NB_MAKE_OPAQUE(Eigen::SparseQRMatrixQReturnType) +NB_MAKE_OPAQUE(Eigen::SparseQRMatrixQTransposeReturnType) +NB_MAKE_OPAQUE(Eigen::SparseLUMatrixLReturnType) +NB_MAKE_OPAQUE(Eigen::SparseLUMatrixUReturnType) + NB_MAKE_OPAQUE(Eigen::LLT) NB_MAKE_OPAQUE(Eigen::LDLT) -NB_MAKE_OPAQUE(Eigen::HouseholderQR) -NB_MAKE_OPAQUE(Eigen::FullPivHouseholderQR) +NB_MAKE_OPAQUE(Eigen::FullPivLU) +NB_MAKE_OPAQUE(Eigen::PartialPivLU) NB_MAKE_OPAQUE(Eigen::ColPivHouseholderQR) NB_MAKE_OPAQUE(Eigen::CompleteOrthogonalDecomposition) +NB_MAKE_OPAQUE(Eigen::FullPivHouseholderQR) +NB_MAKE_OPAQUE(Eigen::HouseholderQR) +NB_MAKE_OPAQUE(Eigen::BDCSVD) +NB_MAKE_OPAQUE(Eigen::ComplexEigenSolver) +NB_MAKE_OPAQUE(Eigen::ComplexSchur) NB_MAKE_OPAQUE(Eigen::EigenSolver) - +NB_MAKE_OPAQUE(Eigen::GeneralizedEigenSolver) +NB_MAKE_OPAQUE(Eigen::GeneralizedSelfAdjointEigenSolver) +NB_MAKE_OPAQUE(Eigen::HessenbergDecomposition) +NB_MAKE_OPAQUE(Eigen::RealQZ) +NB_MAKE_OPAQUE(Eigen::RealSchur) +NB_MAKE_OPAQUE(Eigen::SelfAdjointEigenSolver) +NB_MAKE_OPAQUE(Eigen::Tridiagonalization) + +// Utils std::string printEigenVersion(const char* delim = ".") { std::ostringstream oss; oss << EIGEN_WORLD_VERSION << delim << EIGEN_MAJOR_VERSION << delim @@ -29,42 +72,119 @@ std::string printEigenVersion(const char* delim = ".") { return oss.str(); } -void exposeSolvers(nb::module_& m); - +// Module NB_MODULE(nanoeigenpy, m) { + // exposeConstants(m); + exposePermutationMatrix(m, "PermutationMatrix"); - // Decompositions - exposeLLTSolver(m, "LLT"); - exposeLDLTSolver(m, "LDLT"); - exposeHouseholderQRSolver(m, "HouseholderQR"); - exposeFullPivHouseholderQRSolver(m, "FullPivHouseholderQR"); - exposeColPivHouseholderQRSolver(m, "ColPivHouseholderQR"); - exposeCompleteOrthogonalDecompositionSolver( + // + exposeLDLT(m, "LDLT"); + exposeLLT(m, "LLT"); + // + exposeFullPivLU(m, "FullPivLU"); + exposePartialPivLU(m, "PartialPivLU"); + // + exposeColPivHouseholderQR(m, "ColPivHouseholderQR"); + exposeCompleteOrthogonalDecomposition( m, "CompleteOrthogonalDecomposition"); + exposeFullPivHouseholderQR(m, "FullPivHouseholderQR"); + exposeHouseholderQR(m, "HouseholderQR"); + // + exposeBDCSVD(m, "BDCSVD"); + exposeJacobiSVD(m, "ColPivHhJacobiSVD"); + exposeJacobiSVD(m, "FullPivHhJacobiSVD"); + exposeJacobiSVD(m, "HhJacobiSVD"); + exposeJacobiSVD(m, "NoPrecondJacobiSVD"); + // + exposeComplexEigenSolver(m, "ComplexEigenSolver"); + exposeComplexSchur(m, "ComplexSchur"); exposeEigenSolver(m, "EigenSolver"); + exposeGeneralizedEigenSolver(m, "GeneralizedEigenSolver"); + exposeGeneralizedSelfAdjointEigenSolver( + m, "GeneralizedSelfAdjointEigenSolver"); + exposeHessenbergDecomposition(m, "HessenbergDecomposition"); + exposeRealQZ(m, "RealQZ"); + exposeRealSchur(m, "RealSchur"); exposeSelfAdjointEigenSolver(m, "SelfAdjointEigenSolver"); - exposePermutationMatrix(m, "PermutationMatrix"); + exposeTridiagonalization(m, "Tridiagonalization"); - exposeSimplicialLLT(m, "SimplicialLLT"); + // exposeSimplicialLDLT(m, "SimplicialLDLT"); - + exposeSimplicialLLT(m, "SimplicialLLT"); + // + exposeSparseLU(m, "SparseLU"); + // + exposeSparseQR(m, "SparseQR"); #ifdef NANOEIGENPY_HAS_CHOLMOD + // exposeCholmodSimplicialLLT(m, "CholmodSimplicialLLT"); exposeCholmodSimplicialLDLT(m, "CholmodSimplicialLDLT"); exposeCholmodSupernodalLLT(m, "CholmodSupernodalLLT"); #endif #ifdef NANOEIGENPY_HAS_ACCELERATE + // exposeAccelerate(m); #endif - // Geometry + // exposeQuaternion(m, "Quaternion"); exposeAngleAxis(m, "AngleAxis"); + exposeHyperplane(m, "Hyperplane"); + exposeParametrizedLine(m, "ParametrizedLine"); + exposeRotation2D(m, "Rotation2D"); + exposeUniformScaling(m, "UniformScaling"); + exposeTranslation(m, "Translation"); + + // + exposeJacobiRotation(m, "JacobiRotation"); + + // + nb::module_ solvers = + m.def_submodule("solvers", "Iterative linear solvers in Eigen."); + exposeIdentityPreconditioner(solvers, "IdentityPreconditioner"); + exposeDiagonalPreconditioner(solvers, "DiagonalPreconditioner"); +#if EIGEN_VERSION_AT_LEAST(3, 3, 5) + exposeLeastSquareDiagonalPreconditioner( + solvers, "LeastSquareDiagonalPreconditioner"); +#endif - // Preconditioners (and solvers) - nb::module_ solvers = m.def_submodule("solvers", "Solvers in Eigen."); - exposeSolvers(solvers); + using Eigen::Lower; + + using Eigen::BiCGSTAB; + using Eigen::ConjugateGradient; + using Eigen::DiagonalPreconditioner; + using Eigen::IdentityPreconditioner; + using Eigen::LeastSquareDiagonalPreconditioner; + using Eigen::LeastSquaresConjugateGradient; + using Eigen::MINRES; + + using IdentityConjugateGradient = + ConjugateGradient; + using IdentityLeastSquaresConjugateGradient = + LeastSquaresConjugateGradient; + using DiagonalLeastSquaresConjugateGradient = + LeastSquaresConjugateGradient>; + using IdentityBiCGSTAB = BiCGSTAB; + using DiagonalMINRES = MINRES>; + + exposeConjugateGradient>( + solvers, "ConjugateGradient"); + exposeConjugateGradient( + solvers, "IdentityConjugateGradient"); + exposeLeastSquaresConjugateGradient>( + solvers, "LeastSquaresConjugateGradient"); + exposeLeastSquaresConjugateGradient( + solvers, "IdentityLeastSquaresConjugateGradient"); + exposeLeastSquaresConjugateGradient( + solvers, "DiagonalLeastSquaresConjugateGradient"); + exposeMINRES>(solvers, "MINRES"); + exposeMINRES(solvers, "DiagonalMINRES"); + exposeBiCGSTAB>(solvers, "BiCGSTAB"); + exposeBiCGSTAB(solvers, "IdentityBiCGSTAB"); + + exposeIncompleteLUT(solvers, "IncompleteLUT"); + exposeIncompleteCholesky(solvers, "IncompleteCholesky"); // Utils exposeIsApprox(m); diff --git a/src/solvers.cpp b/src/solvers.cpp deleted file mode 100644 index 31f2310..0000000 --- a/src/solvers.cpp +++ /dev/null @@ -1,35 +0,0 @@ -#include "nanoeigenpy/solvers.hpp" - -#include "./internal.h" - -using namespace nanoeigenpy; - -void exposeSolvers(nb::module_& m) { - exposeIdentityPreconditioner(m, "IdentityPreconditioner"); - exposeDiagonalPreconditioner(m, "DiagonalPreconditioner"); -#if EIGEN_VERSION_AT_LEAST(3, 3, 5) - exposeLeastSquareDiagonalPreconditioner( - m, "LeastSquareDiagonalPreconditioner"); -#endif - exposeMINRESSolver(m, "MINRES"); - - // Solvers - using Eigen::ConjugateGradient; - using Eigen::IdentityPreconditioner; - using Eigen::LeastSquareDiagonalPreconditioner; - using Eigen::LeastSquaresConjugateGradient; - using Eigen::Lower; - using Eigen::Upper; - - exposeConjugateGradient>( - m, "ConjugateGradient"); - - exposeLeastSquaresConjugateGradient>>( - m, "LeastSquaresConjugateGradient"); - - using IdentityConjugateGradient = - ConjugateGradient; - exposeConjugateGradient( - m, "IdentityConjugateGradient"); -} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f15738e..18cfe81 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,3 +1,5 @@ +# Copyright 2025 INRIA + # Create a shared nanobind library for testing set(NANOBIND_TESTING_TARGET nanobind-testing) nanobind_build_library(${NANOBIND_TESTING_TARGET} SHARED) @@ -59,14 +61,30 @@ add_tests_cpp_extension(quaternion) set( TEST_NAMES test_eigen_solver - test_geometry + test_complex_eigen_solver + test_generalized_eigen_solver + test_self_adjoint_eigen_solver + test_generalized_self_adjoint_eigen_solver + test_real_schur + test_complex_schur + test_hessenberg_decomposition + test_real_qz + test_tridiagonalization + test_bdcsvd + test_jacobi_svd + test_full_piv_lu + test_partial_piv_lu test_ldlt test_llt - test_iterative_solvers - test_permutation_matrix test_qr - test_self_adjoint_eigen_solver test_simplicial_llt + test_sparse_lu + test_sparse_qr + test_geometry + test_iterative_solvers + test_permutation_matrix + test_incomplete_lut + test_incomplete_cholesky ) if(BUILD_WITH_CHOLMOD_SUPPORT) diff --git a/tests/quaternion.cpp b/tests/quaternion.cpp index 260d732..82fb99f 100644 --- a/tests/quaternion.cpp +++ b/tests/quaternion.cpp @@ -1,4 +1,5 @@ /// Copyright 2025 INRIA + #include #include diff --git a/tests/test_accelerate.py b/tests/test_accelerate.py index 8a5f6bf..d21d46d 100644 --- a/tests/test_accelerate.py +++ b/tests/test_accelerate.py @@ -19,7 +19,7 @@ def test(SolverType: type): X = rng.random((dim, 20)) B = A.dot(X) X_est = llt.solve(B) - # import pdb; pdb.set_trace() + assert nanoeigenpy.is_approx(X, X_est) assert nanoeigenpy.is_approx(A.dot(X_est), B) diff --git a/tests/test_bdcsvd.py b/tests/test_bdcsvd.py new file mode 100644 index 0000000..054104b --- /dev/null +++ b/tests/test_bdcsvd.py @@ -0,0 +1,141 @@ +import nanoeigenpy +import numpy as np +import pytest + +_options = [ + 0, + nanoeigenpy.DecompositionOptions.ComputeThinU.value, + nanoeigenpy.DecompositionOptions.ComputeThinV.value, + nanoeigenpy.DecompositionOptions.ComputeFullU.value, + nanoeigenpy.DecompositionOptions.ComputeFullV.value, + nanoeigenpy.DecompositionOptions.ComputeThinU.value + | nanoeigenpy.DecompositionOptions.ComputeThinV.value, + nanoeigenpy.DecompositionOptions.ComputeFullU.value + | nanoeigenpy.DecompositionOptions.ComputeFullV.value, + nanoeigenpy.DecompositionOptions.ComputeThinU.value + | nanoeigenpy.DecompositionOptions.ComputeFullV.value, + nanoeigenpy.DecompositionOptions.ComputeFullU.value + | nanoeigenpy.DecompositionOptions.ComputeThinV.value, +] + + +@pytest.mark.parametrize("options", _options) +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 = nanoeigenpy.BDCSVD(A, options) + assert bdcsvd.info() == nanoeigenpy.ComputationInfo.Success + + if options & ( + nanoeigenpy.DecompositionOptions.ComputeThinU.value + | nanoeigenpy.DecompositionOptions.ComputeFullU.value + ) and options & ( + nanoeigenpy.DecompositionOptions.ComputeThinV.value + | nanoeigenpy.DecompositionOptions.ComputeFullV.value + ): + X = rng.random((dim, 20)) + B = A @ X + X_est = bdcsvd.solve(B) + assert nanoeigenpy.is_approx(X, X_est) + assert nanoeigenpy.is_approx(A @ X_est, B) + + x = rng.random(dim) + b = A @ x + x_est = bdcsvd.solve(b) + assert nanoeigenpy.is_approx(x, x_est) + assert nanoeigenpy.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 + & ( + nanoeigenpy.DecompositionOptions.ComputeThinU.value + | nanoeigenpy.DecompositionOptions.ComputeFullU.value + ) + ) + expected_compute_v = bool( + options + & ( + nanoeigenpy.DecompositionOptions.ComputeThinV.value + | nanoeigenpy.DecompositionOptions.ComputeFullV.value + ) + ) + assert compute_u == expected_compute_u + assert compute_v == expected_compute_v + + if compute_u: + matrixU = bdcsvd.matrixU() + if options & nanoeigenpy.DecompositionOptions.ComputeFullU.value: + assert matrixU.shape == (dim, dim) + elif options & nanoeigenpy.DecompositionOptions.ComputeThinU.value: + assert matrixU.shape == (dim, dim) + assert nanoeigenpy.is_approx(matrixU.T @ matrixU, np.eye(matrixU.shape[1])) + + if compute_v: + matrixV = bdcsvd.matrixV() + if options & nanoeigenpy.DecompositionOptions.ComputeFullV.value: + assert matrixV.shape == (dim, dim) + elif options & nanoeigenpy.DecompositionOptions.ComputeThinV.value: + assert matrixV.shape == (dim, dim) + assert nanoeigenpy.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 nanoeigenpy.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 = nanoeigenpy.BDCSVD() + decomp2 = nanoeigenpy.BDCSVD() + id1 = decomp1.id() + id2 = decomp2.id() + assert id1 != id2 + assert id1 == decomp1.id() + assert id2 == decomp2.id() + + decomp3 = nanoeigenpy.BDCSVD(dim, dim, options) + decomp4 = nanoeigenpy.BDCSVD(dim, dim, options) + id3 = decomp3.id() + id4 = decomp4.id() + assert id3 != id4 + assert id3 == decomp3.id() + assert id4 == decomp4.id() + + +if __name__ == "__main__": + import sys + + sys.exit(pytest.main(sys.argv)) diff --git a/tests/test_complex_eigen_solver.py b/tests/test_complex_eigen_solver.py new file mode 100644 index 0000000..df6527a --- /dev/null +++ b/tests/test_complex_eigen_solver.py @@ -0,0 +1,32 @@ +import nanoeigenpy +import numpy as np + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +es = nanoeigenpy.ComplexEigenSolver(A) +assert es.info() == nanoeigenpy.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 nanoeigenpy.is_approx(AV.real, VD.real) +assert nanoeigenpy.is_approx(AV.imag, VD.imag) + +trace_A = np.trace(A) +trace_D = np.sum(D) +assert abs(trace_A - trace_D.real) < 1e-10 +assert abs(trace_D.imag) < 1e-10 + +ces5 = nanoeigenpy.ComplexEigenSolver(A) +ces6 = nanoeigenpy.ComplexEigenSolver(A) +id5 = ces5.id() +id6 = ces6.id() +assert id5 != id6 +assert id5 == ces5.id() +assert id6 == ces6.id() diff --git a/tests/test_complex_schur.py b/tests/test_complex_schur.py new file mode 100644 index 0000000..da1acb9 --- /dev/null +++ b/tests/test_complex_schur.py @@ -0,0 +1,49 @@ +import nanoeigenpy +import numpy as np + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +cs = nanoeigenpy.ComplexSchur(A) +assert cs.info() == nanoeigenpy.ComputationInfo.Success + +U = cs.matrixU() +T = cs.matrixT() + +A_complex = A.astype(complex) +assert nanoeigenpy.is_approx(A_complex, U @ T @ U.conj().T) +assert nanoeigenpy.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 = nanoeigenpy.ComplexSchur(dim) +cs_triangular.setMaxIterations(1) +result_triangular = cs_triangular.compute(A_triangular) +assert result_triangular.info() == nanoeigenpy.ComputationInfo.Success + +T_triangular = cs_triangular.matrixT() +U_triangular = cs_triangular.matrixU() + +A_triangular_complex = A_triangular.astype(complex) +assert nanoeigenpy.is_approx(T_triangular, A_triangular_complex) +assert nanoeigenpy.is_approx(U_triangular, np.eye(dim, dtype=complex)) + +hess = nanoeigenpy.HessenbergDecomposition(A) +H = hess.matrixH() +Q_hess = hess.matrixQ() + +cs_from_hess = nanoeigenpy.ComplexSchur(dim) +result_from_hess = cs_from_hess.computeFromHessenberg(H, Q_hess, True) +assert result_from_hess.info() == nanoeigenpy.ComputationInfo.Success + +T_from_hess = cs_from_hess.matrixT() +U_from_hess = cs_from_hess.matrixU() + +A_complex = A.astype(complex) +assert nanoeigenpy.is_approx( + A_complex, U_from_hess @ T_from_hess @ U_from_hess.conj().T +) diff --git a/tests/test_eigen_solver.py b/tests/test_eigen_solver.py index cff0f9b..27a9594 100644 --- a/tests/test_eigen_solver.py +++ b/tests/test_eigen_solver.py @@ -1,30 +1,22 @@ import nanoeigenpy import numpy as np -dim = 4 -seed = 1 -rng = np.random.default_rng(seed) +dim = 100 +rng = np.random.default_rng() A = rng.random((dim, dim)) -print("A :") -print(A) -# Tests init es = nanoeigenpy.EigenSolver() es = nanoeigenpy.EigenSolver(dim) es = nanoeigenpy.EigenSolver(A) assert es.info() == nanoeigenpy.ComputationInfo.Success -# Test eigenvectors -# Test eigenvalues V = es.eigenvectors() D = es.eigenvalues() assert nanoeigenpy.is_approx(A.dot(V).real, V.dot(np.diag(D)).real) assert nanoeigenpy.is_approx(A.dot(V).imag, V.dot(np.diag(D)).imag) -# Test nb::init<>() -# Test id es1 = nanoeigenpy.EigenSolver() es2 = nanoeigenpy.EigenSolver() @@ -35,8 +27,6 @@ assert id1 == es1.id() assert id2 == es2.id() -# Test nb::init() -# Test id dim_constructor = 3 es3 = nanoeigenpy.EigenSolver(dim_constructor) diff --git a/tests/test_full_piv_lu.py b/tests/test_full_piv_lu.py new file mode 100644 index 0000000..94e5dcc --- /dev/null +++ b/tests/test_full_piv_lu.py @@ -0,0 +1,109 @@ +import nanoeigenpy +import numpy as np + +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 = nanoeigenpy.FullPivLU(A) + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = fullpivlu.solve(B) +assert nanoeigenpy.is_approx(X, X_est) +assert nanoeigenpy.is_approx(A.dot(X_est), B) + +x = rng.random(dim) +b = A.dot(x) +x_est = fullpivlu.solve(b) +assert nanoeigenpy.is_approx(x, x_est) +assert nanoeigenpy.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 nanoeigenpy.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 nanoeigenpy.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() +image = fullpivlu.image(A) +assert kernel.shape[1] == 1 +assert nanoeigenpy.is_approx(A @ kernel, np.zeros((dim, 1))) +assert image.shape[1] == rank + +inverse = fullpivlu.inverse() +assert nanoeigenpy.is_approx(A @ inverse, np.eye(dim)) +assert nanoeigenpy.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 nanoeigenpy.is_approx(P @ P_inv, np.eye(dim)) +assert nanoeigenpy.is_approx(Q @ Q_inv, np.eye(dim)) +assert nanoeigenpy.is_approx(P_inv @ P, np.eye(dim)) +assert nanoeigenpy.is_approx(Q_inv @ Q, np.eye(dim)) + +rows_rect = 4 +cols_rect = 6 +A_rect = rng.random((rows_rect, cols_rect)) +fullpivlu_rect = nanoeigenpy.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 = nanoeigenpy.FullPivLU() +decomp2 = nanoeigenpy.FullPivLU() +id1 = decomp1.id() +id2 = decomp2.id() +assert id1 != id2 +assert id1 == decomp1.id() +assert id2 == decomp2.id() + +decomp3 = nanoeigenpy.FullPivLU(dim, dim) +decomp4 = nanoeigenpy.FullPivLU(dim, dim) +id3 = decomp3.id() +id4 = decomp4.id() +assert id3 != id4 +assert id3 == decomp3.id() +assert id4 == decomp4.id() diff --git a/tests/test_generalized_eigen_solver.py b/tests/test_generalized_eigen_solver.py new file mode 100644 index 0000000..4be542c --- /dev/null +++ b/tests/test_generalized_eigen_solver.py @@ -0,0 +1,40 @@ +import nanoeigenpy +import numpy as np + +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 = nanoeigenpy.GeneralizedEigenSolver(A, B) +assert ges_matrices.info() == nanoeigenpy.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 nanoeigenpy.is_approx(Av.real, lambda_Bv.real, 1e-6) + assert nanoeigenpy.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 nanoeigenpy.is_approx(alpha_Bv.real, beta_Av.real, 1e-6) + assert nanoeigenpy.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/tests/test_generalized_self_adjoint_eigen_solver.py b/tests/test_generalized_self_adjoint_eigen_solver.py new file mode 100644 index 0000000..7ec5f3c --- /dev/null +++ b/tests/test_generalized_self_adjoint_eigen_solver.py @@ -0,0 +1,85 @@ +import nanoeigenpy +import numpy as np +import pytest + +_options = [ + nanoeigenpy.DecompositionOptions.ComputeEigenvectors.value + | nanoeigenpy.DecompositionOptions.Ax_lBx.value, + nanoeigenpy.DecompositionOptions.EigenvaluesOnly.value + | nanoeigenpy.DecompositionOptions.Ax_lBx.value, + nanoeigenpy.DecompositionOptions.ComputeEigenvectors.value + | nanoeigenpy.DecompositionOptions.ABx_lx.value, + nanoeigenpy.DecompositionOptions.EigenvaluesOnly.value + | nanoeigenpy.DecompositionOptions.ABx_lx.value, + nanoeigenpy.DecompositionOptions.ComputeEigenvectors.value + | nanoeigenpy.DecompositionOptions.BAx_lx.value, + nanoeigenpy.DecompositionOptions.EigenvaluesOnly.value + | nanoeigenpy.DecompositionOptions.BAx_lx.value, +] + + +@pytest.mark.parametrize("options", _options) +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 = nanoeigenpy.GeneralizedSelfAdjointEigenSolver(A, B, options) + assert gsaes.info() == nanoeigenpy.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 & nanoeigenpy.DecompositionOptions.ComputeEigenvectors.value + ) + + if compute_eigenvectors: + V = gsaes.eigenvectors() + assert V.shape == (dim, dim) + + if options & nanoeigenpy.DecompositionOptions.Ax_lBx.value: + for i in range(dim): + v = V[:, i] + lam = D[i] + Av = A @ v + lam_Bv = lam * (B @ v) + assert nanoeigenpy.is_approx(Av, lam_Bv, 1e-6) + + VT_B_V = V.T @ B @ V + assert nanoeigenpy.is_approx(VT_B_V, np.eye(dim), 1e-6) + + elif options & nanoeigenpy.DecompositionOptions.ABx_lx.value: + AB = A @ B + for i in range(dim): + v = V[:, i] + lam = D[i] + ABv = AB @ v + lam_v = lam * v + assert nanoeigenpy.is_approx(ABv, lam_v, 1e-6) + + elif options & nanoeigenpy.DecompositionOptions.BAx_lx.value: + BA = B @ A + for i in range(dim): + v = V[:, i] + lam = D[i] + BAv = BA @ v + lam_v = lam * v + assert nanoeigenpy.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 + + +if __name__ == "__main__": + import sys + + sys.exit(pytest.main(sys.argv)) diff --git a/tests/test_geometry.py b/tests/test_geometry.py index 44c6c86..1fe48ab 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -3,32 +3,30 @@ import nanoeigenpy import quaternion + verbose = True def isapprox(a, b, epsilon=1e-6): if issubclass(a.__class__, np.ndarray) and issubclass(b.__class__, np.ndarray): - return np.allclose(a, b, epsilon) + return nanoeigenpy.is_approx(a, b, epsilon) else: return abs(a - b) < epsilon # --- Quaternion --------------------------------------------------------------- -# Coefficient initialisation verbose and print("[Quaternion] Coefficient initialisation") q = nanoeigenpy.Quaternion(1, 2, 3, 4) q.normalize() assert isapprox(np.linalg.norm(q.coeffs()), q.norm()) assert isapprox(np.linalg.norm(q.coeffs()), 1) -# Coefficient-vector initialisation verbose and print("[Quaternion] Coefficient-vector initialisation") v = np.array([0.5, -0.5, 0.5, 0.5]) for k in range(10000): qv = nanoeigenpy.Quaternion(v) assert isapprox(qv.coeffs(), v) -# Angle axis initialisation verbose and print("[Quaternion] AngleAxis initialisation") r = nanoeigenpy.AngleAxis(q) q2 = nanoeigenpy.Quaternion(r) @@ -42,7 +40,6 @@ def isapprox(a, b, epsilon=1e-6): assert isapprox(Rq.dot(Rq.T), np.eye(3)) assert isapprox(Rr, Rq) -# Rotation matrix initialisation verbose and print("[Quaternion] Rotation Matrix initialisation") qR = nanoeigenpy.Quaternion(Rr) assert q.isApprox(qR) @@ -56,7 +53,6 @@ def isapprox(a, b, epsilon=1e-6): if verbose: print("As expected, caught exception: ", e) -# Test struct owning quaternion, constructed from quaternion x = quaternion.X(q) assert x.a == q @@ -79,3 +75,631 @@ def isapprox(a, b, epsilon=1e-6): R = r.matrix() r2 = nanoeigenpy.AngleAxis(np.dot(R, R)) assert isapprox(r2.angle, r.angle * 2) + +# --- Hyperplane ------------------------------------------------ +verbose and print("[Hyperplane] Normal and point construction") +n = np.array([1.0, 0.0]) +p = np.array([2.0, 3.0]) +h = nanoeigenpy.Hyperplane(n, p) +assert isapprox(h.normal(), n) +assert isapprox(h.absDistance(p), 0.0) +assert h.dim() == 2 + +verbose and print("[Hyperplane] Normal and distance construction") +d = -np.dot(n, p) +h2 = nanoeigenpy.Hyperplane(n, d) +assert isapprox(h.coeffs(), h2.coeffs()) +assert isapprox(h2.offset(), d) + +verbose and print("[Hyperplane] Through two points") +p1 = np.array([0.0, 0.0]) +p2 = np.array([1.0, 1.0]) +h3 = nanoeigenpy.Hyperplane.Through(p1, p2) +assert isapprox(h3.absDistance(p1), 0.0) +assert isapprox(h3.absDistance(p2), 0.0) +assert isapprox(np.linalg.norm(h3.normal()), 1.0) + +verbose and print("[Hyperplane] Through three points") +p1_3d = np.array([1.0, 0.0, 0.0]) +p2_3d = np.array([0.0, 1.0, 0.0]) +p3_3d = np.array([0.0, 0.0, 1.0]) +h4 = nanoeigenpy.Hyperplane.Through(p1_3d, p2_3d, p3_3d) +assert isapprox(h4.absDistance(p1_3d), 0.0) +assert isapprox(h4.absDistance(p2_3d), 0.0) +assert isapprox(h4.absDistance(p3_3d), 0.0) +assert isapprox(np.linalg.norm(h4.normal()), 1.0) +assert h4.dim() == 3 + +verbose and print("[Hyperplane] Distance calculations") +test_point = np.array([1.0, 0.0]) +signed_dist = h3.signedDistance(test_point) +abs_dist = h3.absDistance(test_point) +assert isapprox(abs_dist, abs(signed_dist)) + +verbose and print("[Hyperplane] Projection") +proj = h3.projection(test_point) +assert isapprox(h3.absDistance(proj), 0.0) + +verbose and print("[Hyperplane] Normalization") +h_copy = nanoeigenpy.Hyperplane(np.array([2.0, 0.0]), np.array([1.0, 0.0])) +h_copy.normalize() +assert isapprox(np.linalg.norm(h_copy.normal()), 1.0) + +verbose and print("[Hyperplane] Line intersection") +h_line1 = nanoeigenpy.Hyperplane(np.array([1.0, 0.0]), 0.0) +h_line2 = nanoeigenpy.Hyperplane(np.array([0.0, 1.0]), 0.0) +intersection = h_line1.intersection(h_line2) +assert isapprox(intersection, np.array([0.0, 0.0])) + +verbose and print("[Hyperplane] isApprox") +h5 = nanoeigenpy.Hyperplane(h) +assert h.isApprox(h5) +assert h.isApprox(h5, 1e-12) + +# --- ParametrizedLine ------------------------------------------------ +verbose and print("[ParametrizedLine] Origin and direction construction") +origin = np.array([1.0, 2.0]) +direction = np.array([1.0, 0.0]) +line = nanoeigenpy.ParametrizedLine(origin, direction) +assert isapprox(line.origin(), origin) +assert isapprox(line.direction(), direction) +assert line.dim() == 2 + +verbose and print("[ParametrizedLine] Default constructor") +line_default = nanoeigenpy.ParametrizedLine() +assert line_default.dim() == 0 + +verbose and print("[ParametrizedLine] Dimension constructor") +line_3d = nanoeigenpy.ParametrizedLine(3) +assert line_3d.dim() == 3 + +verbose and print("[ParametrizedLine] Copy constructor") +line_copy = nanoeigenpy.ParametrizedLine(line) +assert isapprox(line_copy.origin(), line.origin()) +assert isapprox(line_copy.direction(), line.direction()) + +verbose and print("[ParametrizedLine] Construction from 2D hyperplane") +h_2d = nanoeigenpy.Hyperplane(np.array([1.0, 0.0]), 0.0) +line_from_h = nanoeigenpy.ParametrizedLine(h_2d) +assert line_from_h.dim() == 2 +assert isapprox(line_from_h.origin(), np.array([0.0, 0.0])) +assert isapprox(line_from_h.direction(), np.array([0.0, 1.0])) + +verbose and print("[ParametrizedLine] 3D hyperplane should fail") +h_3d = nanoeigenpy.Hyperplane(np.array([1.0, 0.0, 0.0]), 0.0) +try: + line_fail = nanoeigenpy.ParametrizedLine(h_3d) + print("Error, this message should not appear.") +except ValueError as e: + if verbose: + print("As expected, caught exception:", e) + +verbose and print("[ParametrizedLine] Distance calculations") +test_point = np.array([1.0, 0.0]) +line_x_axis = nanoeigenpy.ParametrizedLine(np.array([0.0, 0.0]), np.array([1.0, 0.0])) +distance = line_x_axis.distance(test_point) +squared_distance = line_x_axis.squaredDistance(test_point) +assert isapprox(distance, 0.0) +assert isapprox(squared_distance, 0.0) + +off_line_point = np.array([1.0, 1.0]) +distance_off = line_x_axis.distance(off_line_point) +squared_distance_off = line_x_axis.squaredDistance(off_line_point) +assert isapprox(distance_off, 1.0) +assert isapprox(squared_distance_off, 1.0) +assert isapprox(distance_off * distance_off, squared_distance_off) + +verbose and print("[ParametrizedLine] Projection") +projection = line_x_axis.projection(off_line_point) +assert isapprox(projection, np.array([1.0, 0.0])) +assert isapprox(line_x_axis.distance(projection), 0.0) + +verbose and print("[ParametrizedLine] Intersection with hyperplane") +line_diagonal = nanoeigenpy.ParametrizedLine(np.array([0.0, 0.0]), np.array([1.0, 1.0])) +h_vertical = nanoeigenpy.Hyperplane(np.array([1.0, 0.0]), -1.0) + +intersection_param = line_diagonal.intersectionParameter(h_vertical) +assert isapprox(intersection_param, 1.0) + +intersection_param_old = line_diagonal.intersection(h_vertical) +assert isapprox(intersection_param_old, intersection_param) + +intersection_point = line_diagonal.intersectionPoint(h_vertical) +expected_intersection = np.array([1.0, 1.0]) +assert isapprox(intersection_point, expected_intersection) +assert isapprox(h_vertical.absDistance(intersection_point), 0.0) + +verbose and print("[ParametrizedLine] isApprox") +line1 = nanoeigenpy.ParametrizedLine(np.array([0.0, 0.0]), np.array([1.0, 0.0])) +line2 = nanoeigenpy.ParametrizedLine(np.array([0.0, 0.0]), np.array([1.0, 0.0])) +line3 = nanoeigenpy.ParametrizedLine(np.array([0.0, 0.0]), np.array([0.0, 1.0])) +assert line1.isApprox(line2) +assert line1.isApprox(line2, 1e-12) +assert not line1.isApprox(line3) + +verbose and print("[ParametrizedLine] Parallel lines") +line_parallel1 = nanoeigenpy.ParametrizedLine( + np.array([0.0, 0.0]), np.array([1.0, 0.0]) +) +line_parallel2 = nanoeigenpy.ParametrizedLine( + np.array([0.0, 1.0]), np.array([1.0, 0.0]) +) +assert not line_parallel1.isApprox(line_parallel2) + +test_points = [np.array([i, 0.0]) for i in range(5)] +distances = [line_parallel2.distance(p) for p in test_points] +for d in distances: + assert isapprox(d, 1.0) + +verbose and print("[ParametrizedLine] Through two points") +p0 = np.array([0.0, 0.0]) +p1 = np.array([1.0, 1.0]) +line_through = nanoeigenpy.ParametrizedLine.Through(p0, p1) +direction = line_through.direction() +expected_dir = (p1 - p0) / np.linalg.norm(p1 - p0) +assert isapprox(line_through.origin(), p0) +assert isapprox(np.linalg.norm(direction), 1.0) +assert isapprox(direction, expected_dir) + +# --- Rotation2D ------------------------------------------------ +verbose and print("[Rotation2D] Default constructor") +r_default = nanoeigenpy.Rotation2D() +assert isapprox(r_default.angle, 0.0) + +verbose and print("[Rotation2D] Angle constructor") +angle = np.pi / 4 +r_angle = nanoeigenpy.Rotation2D(angle) +assert isapprox(r_angle.angle, angle) + +verbose and print("[Rotation2D] Copy constructor") +r_copy = nanoeigenpy.Rotation2D(r_angle) +assert isapprox(r_copy.angle, r_angle.angle) +assert r_copy == r_angle + +verbose and print("[Rotation2D] Matrix constructor") +theta = np.pi / 6 +rotation_matrix = np.array( + [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] +) +r_matrix = nanoeigenpy.Rotation2D(rotation_matrix) +assert isapprox(r_matrix.angle, theta) + +verbose and print("[Rotation2D] Angle property") +r_prop = nanoeigenpy.Rotation2D() +new_angle = np.pi / 3 +r_prop.angle = new_angle +assert isapprox(r_prop.angle, new_angle) + +verbose and print("[Rotation2D] smallestPositiveAngle") +r_negative = nanoeigenpy.Rotation2D(-np.pi / 4) +positive_angle = r_negative.smallestPositiveAngle() +assert positive_angle >= 0.0 +assert positive_angle < 2 * np.pi +assert isapprox(positive_angle, 7 * np.pi / 4) + +verbose and print("[Rotation2D] smallestAngle") +r_large = nanoeigenpy.Rotation2D(3 * np.pi) +smallest_angle = r_large.smallestAngle() +assert smallest_angle >= -np.pi +assert smallest_angle <= np.pi +assert isapprox(smallest_angle, np.pi) + +verbose and print("[Rotation2D] Identity") +r_identity = nanoeigenpy.Rotation2D.Identity() +assert isapprox(r_identity.angle, 0.0) + +verbose and print("[Rotation2D] fromRotationMatrix") +r_from_matrix = nanoeigenpy.Rotation2D() +theta2 = np.pi / 2 +matrix2 = np.array( + [[np.cos(theta2), -np.sin(theta2)], [np.sin(theta2), np.cos(theta2)]] +) +r_from_matrix.fromRotationMatrix(matrix2) +assert isapprox(r_from_matrix.angle, theta2) + +verbose and print("[Rotation2D] Rotation composition") +r1 = nanoeigenpy.Rotation2D(np.pi / 4) +r2 = nanoeigenpy.Rotation2D(np.pi / 6) +r_composed = r1 * r2 +expected_angle = np.pi / 4 + np.pi / 6 +assert isapprox(r_composed.angle, expected_angle) + +verbose and print("[Rotation2D] In-place multiplication") +r_inplace = nanoeigenpy.Rotation2D(np.pi / 4) +original_angle = r_inplace.angle +r_inplace *= nanoeigenpy.Rotation2D(np.pi / 6) +assert isapprox(r_inplace.angle, original_angle + np.pi / 6) + +verbose and print("[Rotation2D] Vector rotation") +r_90 = nanoeigenpy.Rotation2D(np.pi / 2) +vec = np.array([1.0, 0.0]) +rotated_vec = r_90 * vec +expected_vec = np.array([0.0, 1.0]) +assert isapprox(rotated_vec, expected_vec) + +vec2 = np.array([1.0, 1.0]) +r_45 = nanoeigenpy.Rotation2D(np.pi / 4) +rotated_vec2 = r_45 * vec2 +expected_vec2 = np.array([0.0, np.sqrt(2)]) +assert isapprox(rotated_vec2, expected_vec2) + +verbose and print("[Rotation2D] Equality operators") +r_eq1 = nanoeigenpy.Rotation2D(np.pi / 3) +r_eq2 = nanoeigenpy.Rotation2D(np.pi / 3) +r_eq3 = nanoeigenpy.Rotation2D(np.pi / 4) + +assert r_eq1 == r_eq2 +assert not (r_eq1 == r_eq3) +assert r_eq1 != r_eq3 +assert not (r_eq1 != r_eq2) + +verbose and print("[Rotation2D] Periodic angles") +r_period1 = nanoeigenpy.Rotation2D(0.0) +r_period2 = nanoeigenpy.Rotation2D(2 * np.pi) +verbose and print("[Rotation2D] isApprox") +r_approx1 = nanoeigenpy.Rotation2D(np.pi / 4) +r_approx2 = nanoeigenpy.Rotation2D(np.pi / 4 + 1e-15) +r_approx3 = nanoeigenpy.Rotation2D(np.pi / 3) + +assert r_approx1.isApprox(r_approx2) +assert r_approx1.isApprox(r_approx2, 1e-12) +assert not r_approx1.isApprox(r_approx3) + +verbose and print("[Rotation2D] slerp") +r_start = nanoeigenpy.Rotation2D(0.0) +r_end = nanoeigenpy.Rotation2D(np.pi / 2) +r_middle = r_start.slerp(0.5, r_end) +assert isapprox(r_middle.angle, np.pi / 4) + +r_slerp_0 = r_start.slerp(0.0, r_end) +r_slerp_1 = r_start.slerp(1.0, r_end) +assert isapprox(r_slerp_0.angle, r_start.angle) +assert isapprox(r_slerp_1.angle, r_end.angle) + +verbose and print("[Rotation2D] Inverse rotation") +try: + r_original = nanoeigenpy.Rotation2D(np.pi / 3) + r_inverse = r_original.inverse() + assert isapprox(r_inverse.angle, -np.pi / 3) + + r_identity_test = r_original * r_inverse + assert isapprox(r_identity_test.angle, 0.0, 1e-12) +except AttributeError: + if verbose: + print("inverse() method not exposed or not available") + +verbose and print("[Rotation2D] Matrix conversion") +try: + r_matrix_test = nanoeigenpy.Rotation2D(np.pi / 6) + matrix = r_matrix_test.matrix() + + assert matrix.shape == (2, 2) + assert isapprox(matrix @ matrix.T, np.eye(2)) + assert isapprox(np.linalg.det(matrix), 1.0) + + expected_matrix = np.array( + [ + [np.cos(np.pi / 6), -np.sin(np.pi / 6)], + [np.sin(np.pi / 6), np.cos(np.pi / 6)], + ] + ) + assert isapprox(matrix, expected_matrix) + +except AttributeError: + if verbose: + print("matrix() method not exposed or not available") + +verbose and print("[Rotation2D] Angle normalization") +r_large_angle = nanoeigenpy.Rotation2D(3 * np.pi) +vec_test = np.array([1.0, 0.0]) +rotated_large = r_large_angle * vec_test +expected_large = np.array([-1.0, 0.0]) +assert isapprox(rotated_large, expected_large) + +# --- UniformScaling ------------------------------------------------ +verbose and print("[UniformScaling] Default constructor") +s_default = nanoeigenpy.UniformScaling() + +verbose and print("[UniformScaling] Factor constructor") +factor = 2.5 +s_factor = nanoeigenpy.UniformScaling(factor) +assert isapprox(s_factor.factor(), factor) + +verbose and print("[UniformScaling] Copy constructor") +s_copy = nanoeigenpy.UniformScaling(s_factor) +assert isapprox(s_copy.factor(), s_factor.factor()) + +verbose and print("[UniformScaling] Factor getter") +s_test = nanoeigenpy.UniformScaling(3.0) +assert isapprox(s_test.factor(), 3.0) + +verbose and print("[UniformScaling] Inverse scaling") +s_original = nanoeigenpy.UniformScaling(4.0) +s_inverse = s_original.inverse() +assert isapprox(s_inverse.factor(), 1.0 / 4.0) + +s_identity_test = s_original * s_inverse +assert isapprox(s_identity_test.factor(), 1.0) + +verbose and print("[UniformScaling] Concatenation of scalings") +s1 = nanoeigenpy.UniformScaling(2.0) +s2 = nanoeigenpy.UniformScaling(3.0) +s_combined = s1 * s2 +assert isapprox(s_combined.factor(), 6.0) + +verbose and print("[UniformScaling] Multiplication with matrix") +s_scale = nanoeigenpy.UniformScaling(2.0) +matrix = np.array([[1.0, 2.0], [3.0, 4.0]]) +scaled_matrix = s_scale * matrix +expected_matrix = matrix * 2.0 +assert isapprox(scaled_matrix, expected_matrix) + +identity = np.eye(3) +s_identity_scale = nanoeigenpy.UniformScaling(5.0) +scaled_identity = s_identity_scale * identity +expected_identity = identity * 5.0 +assert isapprox(scaled_identity, expected_identity) + +verbose and print("[UniformScaling] Multiplication with AngleAxis") +try: + angle_axis = nanoeigenpy.AngleAxis(np.pi / 4, np.array([0.0, 0.0, 1.0])) + s_with_rotation = nanoeigenpy.UniformScaling(2.0) + result_rotation = s_with_rotation * angle_axis + + assert result_rotation.shape == (3, 3) + det = np.linalg.det(result_rotation) + assert isapprox(det, 2.0**3) + +except (AttributeError, NameError): + if verbose: + print("AngleAxis class not available or not exposed") + +verbose and print("[UniformScaling] Multiplication with Quaternion") +try: + quat = nanoeigenpy.Quaternion(1, 0, 0, 0) + s_with_quat = nanoeigenpy.UniformScaling(3.0) + result_quat = s_with_quat * quat + + assert result_quat.shape == (3, 3) + expected_scaled_identity = np.eye(3) * 3.0 + assert isapprox(result_quat, expected_scaled_identity) + +except (AttributeError, NameError): + if verbose: + print("Quaternion class not available or not exposed") + +verbose and print("[UniformScaling] Multiplication with Rotation2D") +try: + rotation_2d = nanoeigenpy.Rotation2D(np.pi / 4) + s_with_rot2d = nanoeigenpy.UniformScaling(2.0) + result_rot2d = s_with_rot2d * rotation_2d + + assert result_rot2d.shape == (2, 2) + det_2d = np.linalg.det(result_rot2d) + assert isapprox(det_2d, 2.0**2) + +except (AttributeError, NameError): + if verbose: + print("Rotation2D class not available or not exposed") + +verbose and print("[UniformScaling] isApprox") +s_approx1 = nanoeigenpy.UniformScaling(2.0) +s_approx2 = nanoeigenpy.UniformScaling(2.0 + 1e-15) +s_approx3 = nanoeigenpy.UniformScaling(3.0) + +assert s_approx1.isApprox(s_approx2) +assert s_approx1.isApprox(s_approx2, 1e-12) +assert not s_approx1.isApprox(s_approx3) + +verbose and print("[UniformScaling] Edge cases") +s_zero = nanoeigenpy.UniformScaling(0.0) +assert isapprox(s_zero.factor(), 0.0) +try: + s_zero_inverse = s_zero.inverse() + if verbose: + print("Zero scaling inverse:", s_zero_inverse.factor()) +except Exception as e: + if verbose: + print("Zero scaling inverse threw exception (expected):", type(e).__name__) + +s_negative = nanoeigenpy.UniformScaling(-2.0) +assert isapprox(s_negative.factor(), -2.0) +s_negative_inverse = s_negative.inverse() +assert isapprox(s_negative_inverse.factor(), -0.5) + +s_small = nanoeigenpy.UniformScaling(1e-10) +s_small_inverse = s_small.inverse() +assert isapprox(s_small_inverse.factor(), 1e10) + +verbose and print("[UniformScaling] Chain operations") +s_chain1 = nanoeigenpy.UniformScaling(2.0) +s_chain2 = nanoeigenpy.UniformScaling(3.0) +s_chain3 = nanoeigenpy.UniformScaling(4.0) + +left_assoc = (s_chain1 * s_chain2) * s_chain3 +right_assoc = s_chain1 * (s_chain2 * s_chain3) +assert isapprox(left_assoc.factor(), right_assoc.factor()) +assert isapprox(left_assoc.factor(), 24.0) + +verbose and print("[UniformScaling] Vector scaling") +s_vector = nanoeigenpy.UniformScaling(2.0) +vector = np.array([[1.0], [2.0], [3.0]]) +identity_3x3 = np.eye(3) +scaled_identity = s_vector * identity_3x3 +scaled_vector = scaled_identity @ vector +expected_vector = vector * 2.0 +assert isapprox(scaled_vector, expected_vector) + +# --- Translation ------------------------------------------------ +verbose and print("[Translation] Default constructor") +t_default = nanoeigenpy.Translation() + +verbose and print("[Translation] 2D constructor with vector") +t_2d = nanoeigenpy.Translation(np.array([1.0, 2.0])) +assert isapprox(t_2d.x, 1.0) +assert isapprox(t_2d.y, 2.0) + +verbose and print("[Translation] 3D constructor with vector") +t_3d = nanoeigenpy.Translation(np.array([1.0, 2.0, 3.0])) +assert isapprox(t_3d.x, 1.0) +assert isapprox(t_3d.y, 2.0) +assert isapprox(t_3d.z, 3.0) + +verbose and print("[Translation] Vector constructor") +vector = np.array([1.5, 2.5, 3.5]) +t_vector = nanoeigenpy.Translation(vector) +assert isapprox(t_vector.x, 1.5) +assert isapprox(t_vector.y, 2.5) +assert isapprox(t_vector.z, 3.5) + +verbose and print("[Translation] Copy constructor") +t_copy = nanoeigenpy.Translation(t_3d) +assert isapprox(t_copy.x, t_3d.x) +assert isapprox(t_copy.y, t_3d.y) +assert isapprox(t_copy.z, t_3d.z) + +verbose and print("[Translation] Property setters") +t_test = nanoeigenpy.Translation(np.array([1.0, 2.0, 3.0])) +t_test.x = 10.0 +t_test.y = 20.0 +t_test.z = 30.0 +assert isapprox(t_test.x, 10.0) +assert isapprox(t_test.y, 20.0) +assert isapprox(t_test.z, 30.0) + +verbose and print("[Translation] Vector and translation getters") +vector_result = t_test.vector() +translation_result = t_test.translation() +assert isapprox(vector_result[0], 10.0) +assert isapprox(translation_result[0], 10.0) + +verbose and print("[Translation] Inverse") +t_original = nanoeigenpy.Translation(np.array([2.0, 3.0, 4.0])) +t_inverse = t_original.inverse() +assert isapprox(t_inverse.x, -2.0) +assert isapprox(t_inverse.y, -3.0) +assert isapprox(t_inverse.z, -4.0) + +verbose and print("[Translation] Concatenation") +t1 = nanoeigenpy.Translation(np.array([1.0, 2.0, 3.0])) +t2 = nanoeigenpy.Translation(np.array([4.0, 5.0, 6.0])) +t_combined = t1 * t2 +assert isapprox(t_combined.x, 5.0) +assert isapprox(t_combined.y, 7.0) +assert isapprox(t_combined.z, 9.0) + +verbose and print("[Translation] isApprox") +t_approx1 = nanoeigenpy.Translation(np.array([1.0, 2.0, 3.0])) +t_approx2 = nanoeigenpy.Translation(np.array([1.0 + 1e-15, 2.0 + 1e-15, 3.0 + 1e-15])) +t_approx3 = nanoeigenpy.Translation(np.array([1.1, 2.1, 3.1])) +assert t_approx1.isApprox(t_approx2) +assert t_approx1.isApprox(t_approx2, 1e-12) +assert not t_approx1.isApprox(t_approx3) + +# --- JacobiRotation --------------------------------------------------------------- +verbose and print("[JacobiRotation] Default constructor") +j = nanoeigenpy.JacobiRotation() +assert hasattr(j, "c") +assert hasattr(j, "s") + +verbose and print("[JacobiRotation] Cosine-sine constructor") +c_val = 0.8 +s_val = 0.6 +j = nanoeigenpy.JacobiRotation(c_val, s_val) +assert isapprox(j.c, c_val) +assert isapprox(j.s, s_val) + +verbose and print("[JacobiRotation] Property access") +j.c = 0.8 +j.s = 0.6 +assert isapprox(j.c, 0.8) +assert isapprox(j.s, 0.6) +norm_squared = j.c**2 + j.s**2 +assert isapprox(norm_squared, 1.0, 1e-12) + +verbose and print("[JacobiRotation] Multiplication operator") +j1 = nanoeigenpy.JacobiRotation(0.8, 0.6) +j2 = nanoeigenpy.JacobiRotation(0.6, 0.8) +j_mult = j1 * j2 +assert hasattr(j_mult, "c") +assert hasattr(j_mult, "s") +norm_mult = j_mult.c**2 + j_mult.s**2 +assert isapprox(norm_mult, 1.0, 1e-12) + +verbose and print("[JacobiRotation] Transpose") +j = nanoeigenpy.JacobiRotation(0.8, 0.6) +j_t = j.transpose() +assert isapprox(j_t.c, j.c) +assert isapprox(j_t.s, -j.s) + +verbose and print("[JacobiRotation] Adjoint") +j = nanoeigenpy.JacobiRotation(0.8, 0.6) +j_adj = j.adjoint() +assert isapprox(j_adj.c, j.c) +assert isapprox(j_adj.s, -j.s) + +verbose and print("[JacobiRotation] Identity property") +j = nanoeigenpy.JacobiRotation(0.8, 0.6) +j_t = j.transpose() +identity = j * j_t +assert isapprox(identity.c, 1.0, 1e-12) +assert isapprox(identity.s, 0.0, 1e-12) + +verbose and print("[JacobiRotation] makeJacobi from scalars") +j = nanoeigenpy.JacobiRotation() +x, z = 4.0, 1.0 +y = 2.0 +result = j.makeJacobi(x, y, z) +assert isinstance(result, bool) +norm_after = j.c**2 + j.s**2 +assert isapprox(norm_after, 1.0, 1e-12) + +verbose and print("[JacobiRotation] makeJacobi from matrix") +M = np.array([[4.0, 2.0, 1.0], [2.0, 3.0, 0.5], [1.0, 0.5, 1.0]]) +j = nanoeigenpy.JacobiRotation() +result = j.makeJacobi(M, 0, 1) +assert isinstance(result, bool) +norm_matrix = j.c**2 + j.s**2 +assert isapprox(norm_matrix, 1.0, 1e-12) + +verbose and print("[JacobiRotation] makeGivens basic") +j = nanoeigenpy.JacobiRotation() +p_val = 3.0 +q_val = 4.0 +j.makeGivens(p_val, q_val) +norm_givens = j.c**2 + j.s**2 +assert isapprox(norm_givens, 1.0, 1e-12) + +verbose and print("[JacobiRotation] makeGivens with r parameter") +j = nanoeigenpy.JacobiRotation() +p_val = 3.0 +q_val = 4.0 +r_container = np.array([0.0]) +j.makeGivens(p_val, q_val, r_container.ctypes.data) +expected_r = np.sqrt(p_val**2 + q_val**2) + +verbose and print("[JacobiRotation] Edge cases") +j_zero = nanoeigenpy.JacobiRotation(1.0, 0.0) +assert isapprox(j_zero.c, 1.0) +assert isapprox(j_zero.s, 0.0) + +j_90 = nanoeigenpy.JacobiRotation(0.0, 1.0) +assert isapprox(j_90.c, 0.0) +assert isapprox(j_90.s, 1.0) + +j = nanoeigenpy.JacobiRotation() +j.makeGivens(5.0, 0.0) +assert isapprox(abs(j.c), 1.0) +assert isapprox(j.s, 0.0) + +j.makeGivens(0.0, 5.0) +assert isapprox(j.c, 0.0) +assert isapprox(abs(j.s), 1.0) + +verbose and print("[JacobiRotation] makeJacobi small off-diagonal") +j = nanoeigenpy.JacobiRotation() +result = j.makeJacobi(1.0, 1e-15, 2.0) +assert isinstance(result, bool) +if not result: + assert isapprox(j.c, 1.0) + assert isapprox(j.s, 0.0) diff --git a/tests/test_hessenberg_decomposition.py b/tests/test_hessenberg_decomposition.py new file mode 100644 index 0000000..9a09802 --- /dev/null +++ b/tests/test_hessenberg_decomposition.py @@ -0,0 +1,70 @@ +import nanoeigenpy +import numpy as np + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) + +hess = nanoeigenpy.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 nanoeigenpy.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 nanoeigenpy.is_approx(QQ_conj, np.eye(dim)) + +A_test = rng.random((dim, dim)) +hess1 = nanoeigenpy.HessenbergDecomposition(dim) +hess1.compute(A_test) +hess2 = nanoeigenpy.HessenbergDecomposition(A_test) + +H1 = hess1.matrixH() +H2 = hess2.matrixH() +Q1 = hess1.matrixQ() +Q2 = hess2.matrixQ() + +assert nanoeigenpy.is_approx(H1, H2) +assert nanoeigenpy.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 = nanoeigenpy.HessenbergDecomposition(dim) +hess_matrix = nanoeigenpy.HessenbergDecomposition(A) + +hess1_id = nanoeigenpy.HessenbergDecomposition(dim) +hess2_id = nanoeigenpy.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 = nanoeigenpy.HessenbergDecomposition(A) +hess4_id = nanoeigenpy.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/tests/test_incomplete_cholesky.py b/tests/test_incomplete_cholesky.py new file mode 100644 index 0000000..01f4f11 --- /dev/null +++ b/tests/test_incomplete_cholesky.py @@ -0,0 +1,71 @@ +import numpy as np +from scipy.sparse import csc_matrix +import nanoeigenpy + + +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 = nanoeigenpy.solvers.IncompleteCholesky(A) +assert ichol.info() == nanoeigenpy.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() == nanoeigenpy.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/tests/test_incomplete_lut.py b/tests/test_incomplete_lut.py new file mode 100644 index 0000000..4405e93 --- /dev/null +++ b/tests/test_incomplete_lut.py @@ -0,0 +1,49 @@ +import numpy as np +from scipy.sparse import csc_matrix +import nanoeigenpy + +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 = nanoeigenpy.solvers.IncompleteLUT(A) +assert ilut.info() == nanoeigenpy.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() == nanoeigenpy.ComputationInfo.Success + +ilut_params = nanoeigenpy.solvers.IncompleteLUT(A, 1e-4, 15) +assert ilut_params.info() == nanoeigenpy.ComputationInfo.Success + +ilut_set = nanoeigenpy.solvers.IncompleteLUT() +ilut_set.setDroptol(1e-3) +ilut_set.setFillfactor(20) +ilut_set.compute(A) +assert ilut_set.info() == nanoeigenpy.ComputationInfo.Success diff --git a/tests/test_iterative_solvers.py b/tests/test_iterative_solvers.py index 1fbde18..35f7c39 100644 --- a/tests/test_iterative_solvers.py +++ b/tests/test_iterative_solvers.py @@ -3,35 +3,32 @@ import pytest dim = 100 -seed = 6 -rng = np.random.default_rng(seed) +rng = np.random.default_rng() MAX_ITER = 8000 -_clazzes = [ +_classes = [ nanoeigenpy.solvers.ConjugateGradient, nanoeigenpy.solvers.IdentityConjugateGradient, nanoeigenpy.solvers.LeastSquaresConjugateGradient, - nanoeigenpy.solvers.MINRES, + nanoeigenpy.solvers.IdentityLeastSquaresConjugateGradient, + nanoeigenpy.solvers.DiagonalLeastSquaresConjugateGradient, ] -@pytest.mark.parametrize("cls", _clazzes) +@pytest.mark.parametrize("cls", _classes) def test_solver(cls): Q = rng.standard_normal((dim, dim)) A = 0.5 * (Q.T + Q) solver = cls(A) solver.setMaxIterations(MAX_ITER) - # Vector rhs - x = rng.random(dim) b = A.dot(x) x_est = solver.solve(b) + assert solver.info() == nanoeigenpy.ComputationInfo.Success assert nanoeigenpy.is_approx(b, A.dot(x_est), 1e-6) - # Matrix rhs - X = rng.random((dim, 20)) B = A.dot(X) X_est = solver.solve(B) @@ -39,16 +36,13 @@ def test_solver(cls): assert nanoeigenpy.is_approx(B, A.dot(X_est), 1e-6) -@pytest.mark.parametrize("cls", _clazzes) +@pytest.mark.parametrize("cls", _classes) def test_solver_with_guess(cls): Q = rng.standard_normal((dim, dim)) A = 0.5 * (Q.T + Q) solver = cls(A) solver.setMaxIterations(MAX_ITER) - # With guess - # Vector rhs - x = rng.random(dim) b = A.dot(x) x_est = solver.solveWithGuess(b, x + 0.01) @@ -57,8 +51,6 @@ def test_solver_with_guess(cls): assert nanoeigenpy.is_approx(x, x_est, 1e-6) assert nanoeigenpy.is_approx(b, A.dot(x_est), 1e-6) - # Matrix rhs - X = rng.random((dim, 20)) B = A.dot(X) X_est = solver.solveWithGuess(B, X + 0.01) diff --git a/tests/test_jacobi_svd.py b/tests/test_jacobi_svd.py new file mode 100644 index 0000000..8d0f066 --- /dev/null +++ b/tests/test_jacobi_svd.py @@ -0,0 +1,119 @@ +import nanoeigenpy +import numpy as np +import pytest + +THIN_U = nanoeigenpy.DecompositionOptions.ComputeThinU.value +THIN_V = nanoeigenpy.DecompositionOptions.ComputeThinV.value +FULL_U = nanoeigenpy.DecompositionOptions.ComputeFullU.value +FULL_V = nanoeigenpy.DecompositionOptions.ComputeFullV.value + +_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 = [ + nanoeigenpy.ColPivHhJacobiSVD, + # nanoeigenpy.FullPivHhJacobiSVD, + # nanoeigenpy.HhJacobiSVD, + # nanoeigenpy.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, options): + if cls == nanoeigenpy.FullPivHhJacobiSVD: + has_thin_u = bool(options & THIN_U) + has_thin_v = bool(options & THIN_V) + + if has_thin_u or has_thin_v: + return False + + return True + + +@pytest.mark.parametrize("cls", _classes) +@pytest.mark.parametrize("options", _options) +def test_jacobi(cls, options): + if not is_valid_combination(cls, options): + pytest.skip(f"Invalid combination: {cls.__name__} with options {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)) + + jacobisvd = cls(A, options) + assert jacobisvd.info() == nanoeigenpy.ComputationInfo.Success + + has_u = options & (THIN_U | FULL_U) + has_v = options & (THIN_V | FULL_V) + + if has_u and has_v: + X = rng.random((dim, 20)) + B = A @ X + X_est = jacobisvd.solve(B) + assert nanoeigenpy.is_approx(X, X_est) + assert nanoeigenpy.is_approx(A @ X_est, B) + + x = rng.random(dim) + b = A @ x + x_est = jacobisvd.solve(b) + assert nanoeigenpy.is_approx(x, x_est) + assert nanoeigenpy.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, options) + + 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 nanoeigenpy.is_approx(matrixU.T @ matrixU, np.eye(matrixU.shape[1])) + + if compute_v: + matrixV = jacobisvd.matrixV() + assert matrixV.shape == (dim, dim) + assert nanoeigenpy.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 nanoeigenpy.is_approx(A, A_reconstructed) + + +if __name__ == "__main__": + import sys + + sys.exit(pytest.main(sys.argv)) diff --git a/tests/test_ldlt.py b/tests/test_ldlt.py index a511a1d..d844863 100644 --- a/tests/test_ldlt.py +++ b/tests/test_ldlt.py @@ -2,33 +2,24 @@ import numpy as np dim = 100 -seed = 1 -rng = np.random.default_rng(seed) +rng = np.random.default_rng() -# Test isNegative A_neg = -np.eye(dim) ldlt_neg = nanoeigenpy.LDLT(A_neg) assert ldlt_neg.isNegative() assert not ldlt_neg.isPositive() -# Test isPositive A_pos = np.eye(dim) ldlt_pos = nanoeigenpy.LDLT(A_pos) assert ldlt_pos.isPositive() assert not ldlt_pos.isNegative() -# Test nb::init() A = rng.random((dim, dim)) A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim)) ldlt = nanoeigenpy.LDLT(A) - -# Test info assert ldlt.info() == nanoeigenpy.ComputationInfo.Success -# Test matrixL -# Test vector D -# Test transpositionsP L = ldlt.matrixL() D = ldlt.vectorD() P = ldlt.transpositionsP() @@ -36,38 +27,31 @@ np.transpose(P).dot(L.dot(np.diag(D).dot(np.transpose(L).dot(P)))), A ) -# Test solve (matrix) X = rng.random((dim, 20)) B = A.dot(X) X_est = ldlt.solve(B) assert nanoeigenpy.is_approx(X, X_est) assert nanoeigenpy.is_approx(A.dot(X_est), B) -# Test solve (vector) x = rng.random(dim) b = A.dot(x) x_est = ldlt.solve(b) assert nanoeigenpy.is_approx(x, x_est) assert nanoeigenpy.is_approx(A.dot(x_est), b) -# Test reconstructedMatrix A_reconstructed = ldlt.reconstructedMatrix() assert nanoeigenpy.is_approx(A_reconstructed, A) -# Test adjoint adjoint = ldlt.adjoint() assert adjoint is ldlt -# Test rcond A_cond = np.eye(dim) ldlt_cond = nanoeigenpy.LDLT(A_cond) estimated_r_cond_num = ldlt_cond.rcond() assert abs(estimated_r_cond_num - 1) <= 1e-9 -# Test compute ldlt_compute = ldlt.compute(A) -# Test matrixLDLT LDLT = ldlt.matrixLDLT() LDLT_lower_without_diag = np.tril(LDLT, k=-1) L_lower_without_diag = np.tril(L, k=-1) @@ -80,7 +64,6 @@ LDLT_diag = np.diagonal(LDLT) assert nanoeigenpy.is_approx(LDLT_diag, D) -# Test rankUpdate sigma = 3 w = np.ones(dim) ldlt.rankUpdate(w, sigma) @@ -90,8 +73,6 @@ A_updated = np.transpose(P).dot(L.dot(np.diag(D).dot(np.transpose(L).dot(P)))) assert nanoeigenpy.is_approx(A_updated, A + sigma * w * np.transpose(w)) -# Test nb::init<>() -# Test id ldlt1 = nanoeigenpy.LDLT() ldlt2 = nanoeigenpy.LDLT() @@ -102,8 +83,6 @@ assert id1 == ldlt1.id() assert id2 == ldlt2.id() -# Test nb::init() -# Test id dim_constructor = 3 ldlt3 = nanoeigenpy.LDLT(dim_constructor) diff --git a/tests/test_llt.py b/tests/test_llt.py index ae55729..f2b953d 100644 --- a/tests/test_llt.py +++ b/tests/test_llt.py @@ -2,42 +2,34 @@ import numpy as np dim = 100 -seed = 1 -rng = np.random.default_rng(seed) +rng = np.random.default_rng() A = rng.random((dim, dim)) A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim)) -# Test nb::init() llt = nanoeigenpy.LLT(A) -# Test info assert llt.info() == nanoeigenpy.ComputationInfo.Success -# Test matrixL L = llt.matrixL() assert nanoeigenpy.is_approx(L.dot(np.transpose(L)), A) -# Test matrixU U = llt.matrixU() LU = L @ U assert nanoeigenpy.is_approx(LU, A) -# Test solve (matrix) X = rng.random((dim, 20)) B = A.dot(X) X_est = llt.solve(B) assert nanoeigenpy.is_approx(X, X_est) assert nanoeigenpy.is_approx(A.dot(X_est), B) -# Test solve (vector) x = rng.random(dim) b = A.dot(x) x_est = llt.solve(b) assert nanoeigenpy.is_approx(x, x_est) assert nanoeigenpy.is_approx(A.dot(x_est), b) -# Test matrixLLT LLT = llt.matrixLLT() LLT_lower = np.tril(LLT) assert nanoeigenpy.is_approx(LLT_lower, L) @@ -46,24 +38,17 @@ LLT_upper = np.triu(LLT, k=1) assert nanoeigenpy.is_approx(A_upper, LLT_upper) -# Test reconstructedMatrix A_reconstructed = llt.reconstructedMatrix() assert nanoeigenpy.is_approx(A_reconstructed, A) -# Test adjoint adjoint = llt.adjoint() assert adjoint is llt -# Test rcond A_cond = np.eye(dim) llt_cond = nanoeigenpy.LLT(A_cond) estimated_r_cond_num = llt_cond.rcond() assert abs(estimated_r_cond_num - 1) <= 1e-9 -# Test compute -# Done implicitly at init - -# Test rankUpdate sigma = 3 w = np.ones(dim) llt.rankUpdate(w, sigma) @@ -72,8 +57,6 @@ LU = L @ U assert nanoeigenpy.is_approx(LU, A + sigma * w * np.transpose(w)) -# Test nb::init<>() -# Test id llt1 = nanoeigenpy.LLT() llt2 = nanoeigenpy.LLT() @@ -84,8 +67,6 @@ assert id1 == llt1.id() assert id2 == llt2.id() -# Test nb::init() -# Test id dim_constructor = 3 llt3 = nanoeigenpy.LLT(dim_constructor) diff --git a/tests/test_partial_piv_lu.py b/tests/test_partial_piv_lu.py new file mode 100644 index 0000000..98554c6 --- /dev/null +++ b/tests/test_partial_piv_lu.py @@ -0,0 +1,75 @@ +import nanoeigenpy +import numpy as np + +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 = nanoeigenpy.PartialPivLU(A) + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = partialpivlu.solve(B) +assert nanoeigenpy.is_approx(X, X_est) +assert nanoeigenpy.is_approx(A.dot(X_est), B) + +x = rng.random(dim) +b = A.dot(x) +x_est = partialpivlu.solve(b) +assert nanoeigenpy.is_approx(x, x_est) +assert nanoeigenpy.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 nanoeigenpy.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 nanoeigenpy.is_approx(P @ A, L @ U) + +inverse = partialpivlu.inverse() +assert nanoeigenpy.is_approx(A @ inverse, np.eye(dim)) +assert nanoeigenpy.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 nanoeigenpy.is_approx(P @ P_inv, np.eye(dim)) +assert nanoeigenpy.is_approx(P_inv @ P, np.eye(dim)) + +decomp1 = nanoeigenpy.PartialPivLU() +decomp2 = nanoeigenpy.PartialPivLU() +id1 = decomp1.id() +id2 = decomp2.id() +assert id1 != id2 +assert id1 == decomp1.id() +assert id2 == decomp2.id() + +decomp3 = nanoeigenpy.PartialPivLU(dim) +decomp4 = nanoeigenpy.PartialPivLU(dim) +id3 = decomp3.id() +id4 = decomp4.id() +assert id3 != id4 +assert id3 == decomp3.id() +assert id4 == decomp4.id() + +decomp5 = nanoeigenpy.PartialPivLU(A) +decomp6 = nanoeigenpy.PartialPivLU(A) +id5 = decomp5.id() +id6 = decomp6.id() +assert id5 != id6 +assert id5 == decomp5.id() +assert id6 == decomp6.id() diff --git a/tests/test_permutation_matrix.py b/tests/test_permutation_matrix.py index 5a179fe..29e76f8 100644 --- a/tests/test_permutation_matrix.py +++ b/tests/test_permutation_matrix.py @@ -1,34 +1,41 @@ import nanoeigenpy import numpy as np -dim = 10 -seed = 1 -rng = np.random.default_rng(seed) +dim = 100 +rng = np.random.default_rng() indices = rng.permutation(dim) -# Tests init perm = nanoeigenpy.PermutationMatrix(dim) perm = nanoeigenpy.PermutationMatrix(indices) -# Test indices est_indices = perm.indices() assert est_indices.all() == indices.all() -# Test applyTranspositionOnTheLeft -# Test applyTranspositionOnThtRight perm_left = perm.applyTranspositionOnTheLeft(0, 1) perm_left_right = perm_left.applyTranspositionOnTheRight(0, 1) assert perm_left_right.indices().all() == perm.indices().all() -# Test setIdentity perm.setIdentity() assert perm.indices().all() == np.arange(dim).all() dim = dim + 1 perm.setIdentity(dim) assert perm.indices().all() == np.arange(dim).all() -# Test nb::init() -# Test id +perm.setIdentity() +dense = perm.toDenseMatrix() +assert dense.all() == np.eye(dim).all() + +perm = nanoeigenpy.PermutationMatrix(np.array([1, 0, 2])) +perm_t = perm.transpose() +dense = perm.toDenseMatrix() +dense_t = perm_t.toDenseMatrix() +assert dense_t.all() == dense.T.all() + +perm_inv = perm.inverse() +result = perm * perm_inv +identity = result.toDenseMatrix() +assert identity.all() == np.eye(3).all() + dim_constructor = 3 perm1 = nanoeigenpy.PermutationMatrix(dim_constructor) @@ -41,8 +48,6 @@ assert id1 == perm1.id() assert id2 == perm2.id() -# Test nb::init() -# Test id es3 = nanoeigenpy.PermutationMatrix(indices) es4 = nanoeigenpy.PermutationMatrix(indices) diff --git a/tests/test_qr.py b/tests/test_qr.py index 1f509e5..fa0f194 100644 --- a/tests/test_qr.py +++ b/tests/test_qr.py @@ -7,7 +7,6 @@ A = rng.random((rows, cols)) -# Test HouseholderQR decomposition householder_qr = nanoeigenpy.HouseholderQR() householder_qr = nanoeigenpy.HouseholderQR(rows, cols) householder_qr = nanoeigenpy.HouseholderQR(A) @@ -24,7 +23,6 @@ y = householder_qr_eye.solve(x) assert (x == y).all() -# Test FullPivHouseholderQR decomposition fullpiv_householder_qr = nanoeigenpy.FullPivHouseholderQR() fullpiv_householder_qr = nanoeigenpy.FullPivHouseholderQR(rows, cols) fullpiv_householder_qr = nanoeigenpy.FullPivHouseholderQR(A) @@ -46,6 +44,7 @@ y = fullpiv_householder_qr.solve(x) assert (x == y).all() +fullpiv_householder_qr.setThreshold() fullpiv_householder_qr.setThreshold(1e-8) assert fullpiv_householder_qr.threshold() == 1e-8 assert nanoeigenpy.is_approx(np.eye(rows, rows), fullpiv_householder_qr.inverse()) @@ -54,7 +53,6 @@ assert fullpiv_householder_qr.nonzeroPivots() == rows assert fullpiv_householder_qr.dimensionOfKernel() == 0 -# Test ColPivHouseholderQR decomposition colpiv_householder_qr = nanoeigenpy.ColPivHouseholderQR(A) assert colpiv_householder_qr.info() == nanoeigenpy.ComputationInfo.Success @@ -67,6 +65,7 @@ assert (X == Y).all() assert colpiv_householder_qr.rank() == rows +colpiv_householder_qr.setThreshold() colpiv_householder_qr.setThreshold(1e-8) assert colpiv_householder_qr.threshold() == 1e-8 assert nanoeigenpy.is_approx(np.eye(rows, rows), colpiv_householder_qr.inverse()) @@ -75,7 +74,6 @@ assert colpiv_householder_qr.nonzeroPivots() == rows assert colpiv_householder_qr.dimensionOfKernel() == 0 -# Test CompleteOrthogonalDecomposition cod = nanoeigenpy.CompleteOrthogonalDecomposition(A) assert cod.info() == nanoeigenpy.ComputationInfo.Success @@ -92,6 +90,7 @@ y = cod.solve(x) assert (x == y).all() +cod.setThreshold() cod.setThreshold(1e-8) assert cod.threshold() == 1e-8 assert nanoeigenpy.is_approx(np.eye(rows, rows), cod.pseudoInverse()) diff --git a/tests/test_real_qz.py b/tests/test_real_qz.py new file mode 100644 index 0000000..088708d --- /dev/null +++ b/tests/test_real_qz.py @@ -0,0 +1,37 @@ +import nanoeigenpy +import numpy as np + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +B = rng.random((dim, dim)) + +realqz = nanoeigenpy.RealQZ(A, B) +assert realqz.info() == nanoeigenpy.ComputationInfo.Success + +Q = realqz.matrixQ() +S = realqz.matrixS() +Z = realqz.matrixZ() +T = realqz.matrixT() + +assert nanoeigenpy.is_approx(A, Q @ S @ Z) +assert nanoeigenpy.is_approx(B, Q @ T @ Z) + +assert nanoeigenpy.is_approx(Q @ Q.T, np.eye(dim)) +assert nanoeigenpy.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 = nanoeigenpy.RealQZ(A, B) +realqz4_id = nanoeigenpy.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/tests/test_real_schur.py b/tests/test_real_schur.py new file mode 100644 index 0000000..399675a --- /dev/null +++ b/tests/test_real_schur.py @@ -0,0 +1,48 @@ +import nanoeigenpy +import numpy as np + + +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 = nanoeigenpy.RealSchur(A) +assert rs.info() == nanoeigenpy.ComputationInfo.Success + +U = rs.matrixU() +T = rs.matrixT() + +assert nanoeigenpy.is_approx(A, U @ T @ U.T) +assert nanoeigenpy.is_approx(U @ U.T, np.eye(dim)) + +verify_is_quasi_triangular(T) + +hess = nanoeigenpy.HessenbergDecomposition(A) +H = hess.matrixH() +Q_hess = hess.matrixQ() + +rs_from_hess = nanoeigenpy.RealSchur(dim) +result_from_hess = rs_from_hess.computeFromHessenberg(H, Q_hess, True) +assert result_from_hess.info() == nanoeigenpy.ComputationInfo.Success + +T_from_hess = rs_from_hess.matrixT() +U_from_hess = rs_from_hess.matrixU() + +assert nanoeigenpy.is_approx(A, U_from_hess @ T_from_hess @ U_from_hess.T) diff --git a/tests/test_self_adjoint_eigen_solver.py b/tests/test_self_adjoint_eigen_solver.py index 4a7ee1a..d53ed89 100644 --- a/tests/test_self_adjoint_eigen_solver.py +++ b/tests/test_self_adjoint_eigen_solver.py @@ -1,9 +1,8 @@ import nanoeigenpy import numpy as np -dim = 5 -seed = 1 -rng = np.random.default_rng(seed) +dim = 100 +rng = np.random.default_rng() A = rng.random((dim, dim)) A = (A + A.T) * 0.5 @@ -18,4 +17,4 @@ AdotV = A @ V VdotD = V @ np.diag(D) -assert nanoeigenpy.is_approx(AdotV, VdotD, 1e-10) +assert nanoeigenpy.is_approx(AdotV, VdotD, 1e-6) diff --git a/tests/test_simplicial_llt.py b/tests/test_simplicial_llt.py index c3878fe..12e8d46 100644 --- a/tests/test_simplicial_llt.py +++ b/tests/test_simplicial_llt.py @@ -2,9 +2,8 @@ import numpy as np import scipy.sparse as spa -dim = 10 -np.set_printoptions(precision=3, linewidth=200) -rng = np.random.default_rng(30) +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 @@ -38,7 +37,6 @@ X_sparse = spa.random(dim, 10, random_state=rng) B_sparse = A.dot(X_sparse) B_sparse: spa.csc_matrix = B_sparse.tocsc(True) -# super important if not B_sparse.has_sorted_indices: B_sparse.sort_indices() diff --git a/tests/test_sparse_lu.py b/tests/test_sparse_lu.py new file mode 100644 index 0000000..859b1f7 --- /dev/null +++ b/tests/test_sparse_lu.py @@ -0,0 +1,63 @@ +import nanoeigenpy +import numpy as np +import scipy.sparse as spa + +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 = nanoeigenpy.SparseLU(A) + +assert splu.info() == nanoeigenpy.ComputationInfo.Success + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = splu.solve(B) +assert isinstance(X_est, np.ndarray) +assert nanoeigenpy.is_approx(X, X_est) +assert nanoeigenpy.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 nanoeigenpy.is_approx(X_est.toarray(), X_sparse.toarray()) +assert nanoeigenpy.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 nanoeigenpy.is_approx(x_reconstructed, x_true, 1e-6) diff --git a/tests/test_sparse_qr.py b/tests/test_sparse_qr.py new file mode 100644 index 0000000..7a825f8 --- /dev/null +++ b/tests/test_sparse_qr.py @@ -0,0 +1,73 @@ +import nanoeigenpy +import numpy as np +import scipy.sparse as spa + +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 = nanoeigenpy.SparseQR(A) + +assert spqr.info() == nanoeigenpy.ComputationInfo.Success + +X = rng.random((dim, 20)) +B = A.dot(X) +X_est = spqr.solve(B) +assert isinstance(X_est, np.ndarray) +assert nanoeigenpy.is_approx(X, X_est) +assert nanoeigenpy.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 nanoeigenpy.is_approx(X_est.toarray(), X_sparse.toarray()) +assert nanoeigenpy.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 nanoeigenpy.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 nanoeigenpy.is_approx(QtAP, R_dense) diff --git a/tests/test_tridiagonalization.py b/tests/test_tridiagonalization.py new file mode 100644 index 0000000..7b8a2f7 --- /dev/null +++ b/tests/test_tridiagonalization.py @@ -0,0 +1,103 @@ +import nanoeigenpy +import numpy as np + +dim = 100 +rng = np.random.default_rng() +A = rng.random((dim, dim)) +A = (A + A.T) * 0.5 + +tri = nanoeigenpy.Tridiagonalization(A) + +Q = tri.matrixQ() +T = tri.matrixT() + +assert nanoeigenpy.is_approx(A, Q @ T @ Q.T) +assert nanoeigenpy.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 nanoeigenpy.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 = nanoeigenpy.Tridiagonalization(dim) +tri1.compute(A_test) +tri2 = nanoeigenpy.Tridiagonalization(A_test) + +Q1 = tri1.matrixQ() +T1 = tri1.matrixT() +Q2 = tri2.matrixQ() +T2 = tri2.matrixT() + +assert nanoeigenpy.is_approx(Q1, Q2) +assert nanoeigenpy.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 = nanoeigenpy.Tridiagonalization(A_diag) +Q_diag = tri_diag.matrixQ() +T_diag = tri_diag.matrixT() + +assert nanoeigenpy.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 = nanoeigenpy.Tridiagonalization(A_tridiag) +Q_tridiag = tri_tridiag.matrixQ() +T_tridiag = tri_tridiag.matrixT() + +assert nanoeigenpy.is_approx(A_tridiag, Q_tridiag @ T_tridiag @ Q_tridiag.T) + + +tri1_id = nanoeigenpy.Tridiagonalization(dim) +tri2_id = nanoeigenpy.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 = nanoeigenpy.Tridiagonalization(A) +tri4_id = nanoeigenpy.Tridiagonalization(A) +id3 = tri3_id.id() +id4 = tri4_id.id() +assert id3 != id4 +assert id3 == tri3_id.id() +assert id4 == tri4_id.id()