diff --git a/cpp/modmesh/linalg/CMakeLists.txt b/cpp/modmesh/linalg/CMakeLists.txt index 00f19c332..2e26dddce 100644 --- a/cpp/modmesh/linalg/CMakeLists.txt +++ b/cpp/modmesh/linalg/CMakeLists.txt @@ -13,6 +13,7 @@ set(MODMESH_LINALG_HEADERS CACHE FILEPATH "" FORCE) set(MODMESH_LINALG_SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/EigenSystem.cpp CACHE FILEPATH "" FORCE) set(MODMESH_LINALG_PYMODHEADERS diff --git a/cpp/modmesh/linalg/EigenSystem.cpp b/cpp/modmesh/linalg/EigenSystem.cpp new file mode 100644 index 000000000..d7174221e --- /dev/null +++ b/cpp/modmesh/linalg/EigenSystem.cpp @@ -0,0 +1,139 @@ +/* + * Copyright (c) 2026, Yung-Yu Chen + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * - Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * - Neither the name of the copyright holder nor the names of its contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include +#include +#include + +// EigenSystem.hpp requires a vendor LAPACK; only compile the surrogate's +// out-of-line bodies when one is available. Without it this is an empty TU. +#ifdef MM_HAS_VENDOR_LAPACK + +#include + +namespace modmesh +{ + +// Reinterpret the plex's type-erased instance as the typed SimpleArray it +// actually holds; the caller must have verified data_type() == T. +template +static SimpleArray const & typed_array(SimpleArrayPlex const & plex) +{ + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast) + return *reinterpret_cast const *>(plex.instance_ptr()); +} + +EigenSystemPlex::EigenSystemPlex(SimpleArrayPlex const & matrix, bool do_vl, bool do_vr) + : m_matrix(matrix) + , m_data_type(matrix.data_type()) +{ + switch (matrix.data_type()) + { + case DataType::Float32: + m_solver = EigenSystem::construct(typed_array(matrix), do_vl, do_vr); + break; + case DataType::Float64: + m_solver = EigenSystem::construct(typed_array(matrix), do_vl, do_vr); + break; + case DataType::Complex64: + m_solver = EigenSystem>::construct(typed_array>(matrix), do_vl, do_vr); + break; + case DataType::Complex128: + m_solver = EigenSystem>::construct(typed_array>(matrix), do_vl, do_vr); + break; + default: + throw std::invalid_argument( + "EigenSystemPlex: SimpleArray data type must be float32, float64, complex64, or complex128"); + } +} + +void EigenSystemPlex::run() +{ + std::visit([](auto const & solver) + { solver->run(); }, + m_solver); +} + +SimpleArrayPlex EigenSystemPlex::wr() const +{ + return std::visit( + [](auto const & solver) + { return SimpleArrayPlex(solver->wr()); }, + m_solver); +} + +SimpleArrayPlex EigenSystemPlex::wi() const +{ + return std::visit( + [](auto const & solver) + { return SimpleArrayPlex(solver->wi()); }, + m_solver); +} + +SimpleArrayPlex EigenSystemPlex::vl(bool suppress_exception) const +{ + return std::visit( + [suppress_exception](auto const & solver) + { return SimpleArrayPlex(solver->vl(suppress_exception)); }, + m_solver); +} + +SimpleArrayPlex EigenSystemPlex::vr(bool suppress_exception) const +{ + return std::visit( + [suppress_exception](auto const & solver) + { return SimpleArrayPlex(solver->vr(suppress_exception)); }, + m_solver); +} + +bool EigenSystemPlex::do_vl() const +{ + return std::visit([](auto const & solver) + { return solver->do_vl(); }, + m_solver); +} + +bool EigenSystemPlex::do_vr() const +{ + return std::visit([](auto const & solver) + { return solver->do_vr(); }, + m_solver); +} + +bool EigenSystemPlex::done() const +{ + return std::visit([](auto const & solver) + { return solver->done(); }, + m_solver); +} + +} /* end namespace modmesh */ + +#endif /* MM_HAS_VENDOR_LAPACK */ + +// vim: set ff=unix fenc=utf8 et sw=4 ts=4 sts=4: diff --git a/cpp/modmesh/linalg/EigenSystem.hpp b/cpp/modmesh/linalg/EigenSystem.hpp index 6252d949d..467b7301e 100644 --- a/cpp/modmesh/linalg/EigenSystem.hpp +++ b/cpp/modmesh/linalg/EigenSystem.hpp @@ -29,8 +29,9 @@ */ /** - * Eigenvalue and eigenvector computation for general (non-symmetric) real - * matrices using LAPACK DGEEV. + * Eigenvalue and eigenvector computation for general (non-symmetric) matrices + * using LAPACK *GEEV: SGEEV (float), DGEEV (double), CGEEV (Complex), + * and ZGEEV (Complex). * * The LAPACK backend is selected by the build system via MM_HAS_VENDOR_LAPACK: * Apple's vecLib (Accelerate framework) on macOS and OpenBLAS on Linux. @@ -43,8 +44,10 @@ #include #include #include +#include #include #include +#include #include #include @@ -53,21 +56,42 @@ namespace modmesh { /** - * Eigenvalue solver for a real general matrix using LAPACK DGEEV. + * Eigenvalue solver for a general matrix using LAPACK *GEEV. * - * Construction validates the input shape and prepares column-major workspace - * buffers. Call run() to invoke DGEEV to calculate eigenvalues and - * eigenvectors. + * The element type T is one of float, double, Complex, or + * Complex. Construction validates the input shape and prepares + * column-major workspace buffers. Call run() to invoke the matching *GEEV + * routine to calculate eigenvalues and eigenvectors. */ +template class EigenSystem + : public std::enable_shared_from_this> { + static_assert( + is_real_v || is_complex_v, + "EigenSystem requires T to be a real or complex number type"); + +private: + + struct ctor_passkey + { + }; + public: - using value_type = double; - using array_type = SimpleArray; + using value_type = T; + using array_type = SimpleArray; + using real_type = typename detail::select_real_t::type; + using real_array_type = SimpleArray; + static constexpr bool is_complex = is_complex_v; - explicit EigenSystem(array_type const & matrix, bool do_vl = true, bool do_vr = true); + static std::shared_ptr construct(array_type const & matrix, bool do_vl = true, bool do_vr = true) + { + return std::make_shared(matrix, do_vl, do_vr, ctor_passkey()); + } + + explicit EigenSystem(array_type const & matrix, bool do_vl, bool do_vr, ctor_passkey const &); EigenSystem() = delete; EigenSystem(EigenSystem const &) = delete; @@ -79,8 +103,8 @@ class EigenSystem void run(); array_type const & matrix() const { return m_matrix; } - array_type const & wr() const { return m_wr; } - array_type const & wi() const { return m_wi; } + real_array_type const & wr() const { return m_wr; } + real_array_type const & wi() const { return m_wi; } array_type const & vl(bool suppress_exception = false) const; array_type const & vr(bool suppress_exception = false) const; bool do_vl() const { return m_do_vl; } @@ -91,19 +115,20 @@ class EigenSystem static std::string format_shape(array_type const & arr); - SimpleArray const & m_matrix; - SimpleArray m_colmajor; - SimpleArray m_wr; - SimpleArray m_wi; - SimpleArray m_vl; - SimpleArray m_vr; + array_type const & m_matrix; + array_type m_colmajor; + real_array_type m_wr; + real_array_type m_wi; + array_type m_vl; + array_type m_vr; bool const m_do_vl; bool const m_do_vr; bool m_done = false; }; /* end class EigenSystem */ -inline EigenSystem::EigenSystem(array_type const & matrix, bool do_vl, bool do_vr) +template +EigenSystem::EigenSystem(array_type const & matrix, bool do_vl, bool do_vr, ctor_passkey const &) : m_matrix(matrix) // Stage matrix into a column-major workspace for LAPACK. , m_colmajor(matrix.to_column_major()) @@ -132,9 +157,10 @@ inline EigenSystem::EigenSystem(array_type const & matrix, bool do_vl, bool do_v } /** - * @brief Run DGEEV on the prepared workspace. + * @brief Run the matching *GEEV routine on the prepared workspace. */ -inline void EigenSystem::run() +template +void EigenSystem::run() { auto const n = static_cast(m_matrix.shape(0)); if (n == 0) @@ -144,88 +170,94 @@ inline void EigenSystem::run() } /* - * Apple Accelerate DGEEV reference: - * https://developer.apple.com/documentation/accelerate/dgeev_(_:_:_:_:_:_:_:_:_:_:_:_:_:_:) - * LAPACK DGEEV API reference: + * LAPACK *GEEV API reference (real DGEEV shown; the complex CGEEV/ZGEEV + * return a single complex eigenvalue array w and take a real rwork): * https://www.netlib.org/lapack/explore-html/d4/d68/group__geev_ga7d8afe93d23c5862e238626905ee145e.html * https://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f_source.html */ char const jobvl = m_do_vl ? 'V' : 'N'; char const jobvr = m_do_vr ? 'V' : 'N'; lapack_int_t const lda = n; - // DGEEV requires LDVL/LDVR >= 1 and a valid (non-null) pointer even - // when the matrix is unreferenced; route the unused side to a stack - // scratch slot. + // *GEEV requires LDVL/LDVR >= 1 and a valid (non-null) pointer even when + // the matrix is unreferenced; route the unused side to a stack scratch + // slot. lapack_int_t const ldvl = m_do_vl ? n : 1; lapack_int_t const ldvr = m_do_vr ? n : 1; lapack_int_t info = 0; - double vl_dummy = 0.0; - double vr_dummy = 0.0; - double * const vl_ptr = m_do_vl ? m_vl.data() : &vl_dummy; - double * const vr_ptr = m_do_vr ? m_vr.data() : &vr_dummy; - - // Phase 1: workspace query. lwork == -1 tells DGEEV to write the - // optimal workspace size into work[0] without performing any work. - double work_query = 0.0; + value_type vl_dummy{}; + value_type vr_dummy{}; + value_type * const vl_ptr = m_do_vl ? m_vl.data() : &vl_dummy; + value_type * const vr_ptr = m_do_vr ? m_vr.data() : &vr_dummy; + + // Phase 1: workspace query. lwork == -1 tells *GEEV to write the optimal + // workspace size into work[0] without performing any work. Phase 2 then + // allocates max(optimal, documented minimum) and runs the factorization. + value_type work_query{}; lapack_int_t lwork = -1; - dgeev_( - &jobvl, - &jobvr, - &n, - m_colmajor.data(), - &lda, - m_wr.data(), - m_wi.data(), - vl_ptr, - &ldvl, - vr_ptr, - &ldvr, - &work_query, - &lwork, - &info); - if (info != 0) + + if constexpr (is_complex) { - throw std::runtime_error(std::format( - "EigenSystem::run: DGEEV workspace query failed with info={}", - static_cast(info))); + // CGEEV/ZGEEV return eigenvalues as a single complex array and need a + // real workspace rwork of length 2*n. Compute into a local complex + // buffer, then split it into the real m_wr/m_wi parts below. + array_type w(static_cast(n)); + real_array_type rwork(static_cast(2 * n)); + detail::lapack_geev( + jobvl, jobvr, n, m_colmajor.data(), lda, w.data(), vl_ptr, ldvl, vr_ptr, ldvr, &work_query, lwork, rwork.data(), &info); + if (info != 0) + { + throw std::runtime_error(std::format( + "EigenSystem::run: GEEV workspace query failed with info={}", + static_cast(info))); + } + + lapack_int_t const lwork_min = 2 * n; + lwork = std::max(static_cast(work_query.real()), lwork_min); + array_type work(static_cast(lwork)); + detail::lapack_geev( + jobvl, jobvr, n, m_colmajor.data(), lda, w.data(), vl_ptr, ldvl, vr_ptr, ldvr, work.data(), lwork, rwork.data(), &info); + for (lapack_int_t i = 0; i < n; ++i) + { + m_wr(i) = w(i).real(); + m_wi(i) = w(i).imag(); + } + } + else + { + detail::lapack_geev( + jobvl, jobvr, n, m_colmajor.data(), lda, m_wr.data(), m_wi.data(), vl_ptr, ldvl, vr_ptr, ldvr, &work_query, lwork, &info); + if (info != 0) + { + throw std::runtime_error(std::format( + "EigenSystem::run: GEEV workspace query failed with info={}", + static_cast(info))); + } + + // 4*n minimum when any eigenvectors requested, else 3*n. + lapack_int_t const lwork_min = (m_do_vl || m_do_vr) ? 4 * n : 3 * n; + lwork = std::max(static_cast(work_query), lwork_min); + array_type work(static_cast(lwork)); + detail::lapack_geev( + jobvl, jobvr, n, m_colmajor.data(), lda, m_wr.data(), m_wi.data(), vl_ptr, ldvl, vr_ptr, ldvr, work.data(), lwork, &info); } - // Phase 2: 4*n minimum when any eigenvectors requested, else 3*n. - lapack_int_t const lwork_min = (m_do_vl || m_do_vr) ? 4 * n : 3 * n; - lwork = std::max(static_cast(work_query), lwork_min); - array_type work(static_cast(lwork)); - dgeev_( - &jobvl, - &jobvr, - &n, - m_colmajor.data(), - &lda, - m_wr.data(), - m_wi.data(), - vl_ptr, - &ldvl, - vr_ptr, - &ldvr, - work.data(), - &lwork, - &info); if (info < 0) { throw std::runtime_error(std::format( - "EigenSystem::run: DGEEV illegal argument at position {}", + "EigenSystem::run: GEEV illegal argument at position {}", -static_cast(info))); } if (info > 0) { throw std::runtime_error(std::format( - "EigenSystem::run: DGEEV failed to converge; {} eigenvalues did not converge", + "EigenSystem::run: GEEV failed to converge; {} eigenvalues did not converge", static_cast(info))); } m_done = true; } /** - * @brief Left eigenvectors computed by DGEEV. + * @brief Left eigenvectors computed by *GEEV. * * @param suppress_exception * Return the empty placeholder instead of throwing when the matrix was @@ -235,7 +267,8 @@ inline void EigenSystem::run() * suppress_exception=true. * @throws std::runtime_error When do_vl=false and suppress_exception=false. */ -inline EigenSystem::array_type const & EigenSystem::vl(bool suppress_exception) const +template +EigenSystem::array_type const & EigenSystem::vl(bool suppress_exception) const { if (!m_do_vl && !suppress_exception) { @@ -247,7 +280,7 @@ inline EigenSystem::array_type const & EigenSystem::vl(bool suppress_exception) } /** - * @brief Right eigenvectors computed by DGEEV. + * @brief Right eigenvectors computed by *GEEV. * * @param suppress_exception * Return the empty placeholder instead of throwing when the matrix was @@ -257,7 +290,8 @@ inline EigenSystem::array_type const & EigenSystem::vl(bool suppress_exception) * suppress_exception=true. * @throws std::runtime_error When do_vr=false and suppress_exception=false. */ -inline EigenSystem::array_type const & EigenSystem::vr(bool suppress_exception) const +template +EigenSystem::array_type const & EigenSystem::vr(bool suppress_exception) const { if (!m_do_vr && !suppress_exception) { @@ -268,7 +302,8 @@ inline EigenSystem::array_type const & EigenSystem::vr(bool suppress_exception) return m_vr; } -inline std::string EigenSystem::format_shape(array_type const & arr) +template +std::string EigenSystem::format_shape(array_type const & arr) { std::string result = "("; for (size_t i = 0; i < arr.ndim(); ++i) @@ -283,6 +318,54 @@ inline std::string EigenSystem::format_shape(array_type const & arr) return result; } +/** + * Type-erased surrogate for EigenSystem. + * + * Accepts a SimpleArrayPlex and dispatches on its runtime element type to the + * matching EigenSystem (float, double, Complex, Complex). + * Results are returned as type-erased SimpleArrayPlex. + */ +class EigenSystemPlex +{ + +public: + + explicit EigenSystemPlex(SimpleArrayPlex const & matrix, bool do_vl = true, bool do_vr = true); + + EigenSystemPlex() = delete; + EigenSystemPlex(EigenSystemPlex const &) = delete; + EigenSystemPlex(EigenSystemPlex &&) = delete; + EigenSystemPlex & operator=(EigenSystemPlex const &) = delete; + EigenSystemPlex & operator=(EigenSystemPlex &&) = delete; + ~EigenSystemPlex() = default; + + void run(); + + SimpleArrayPlex const & matrix() const { return m_matrix; } + DataType data_type() const { return m_data_type; } + + SimpleArrayPlex wr() const; + SimpleArrayPlex wi() const; + SimpleArrayPlex vl(bool suppress_exception = false) const; + SimpleArrayPlex vr(bool suppress_exception = false) const; + bool do_vl() const; + bool do_vr() const; + bool done() const; + +private: + + using solver_variant = std::variant< + std::shared_ptr>, + std::shared_ptr>, + std::shared_ptr>>, + std::shared_ptr>>>; + + SimpleArrayPlex const & m_matrix; + solver_variant m_solver; + DataType m_data_type; + +}; /* end class EigenSystemPlex */ + } /* end namespace modmesh */ // vim: set ff=unix fenc=utf8 et sw=4 ts=4 sts=4: diff --git a/cpp/modmesh/linalg/lapack_compat.hpp b/cpp/modmesh/linalg/lapack_compat.hpp index 51542cd8d..600f2beb5 100644 --- a/cpp/modmesh/linalg/lapack_compat.hpp +++ b/cpp/modmesh/linalg/lapack_compat.hpp @@ -46,6 +46,8 @@ #error "modmesh/linalg/lapack_compat.hpp requires a vendor LAPACK (MM_HAS_VENDOR_LAPACK)." #endif +#include + #ifdef __APPLE__ // Opt in to the modern (non-deprecated) LAPACK signatures provided by // Accelerate. Must be defined before including the Accelerate header. @@ -55,6 +57,26 @@ // NOLINTNEXTLINE(misc-header-include-cycle) #include #else // __APPLE__ +// On Linux + OpenBLAS we declare the Fortran ABI symbols directly. The +// general (non-symmetric) eigensolvers are SGEEV/DGEEV (real) and +// CGEEV/ZGEEV (complex). The complex buffers are passed as void* so this +// header needs no LAPACK complex type; the real workspace rwork is a plain +// float/double array. +extern "C" void sgeev_( + char const * jobvl, + char const * jobvr, + int const * n, + float * a, + int const * lda, + float * wr, + float * wi, + float * vl, + int const * ldvl, + float * vr, + int const * ldvr, + float * work, + int const * lwork, + int * info); extern "C" void dgeev_( char const * jobvl, char const * jobvr, @@ -70,6 +92,36 @@ extern "C" void dgeev_( double * work, int const * lwork, int * info); +extern "C" void cgeev_( + char const * jobvl, + char const * jobvr, + int const * n, + void * a, + int const * lda, + void * w, + void * vl, + int const * ldvl, + void * vr, + int const * ldvr, + void * work, + int const * lwork, + float * rwork, + int * info); +extern "C" void zgeev_( + char const * jobvl, + char const * jobvr, + int const * n, + void * a, + int const * lda, + void * w, + void * vl, + int const * ldvl, + void * vr, + int const * ldvr, + void * work, + int const * lwork, + double * rwork, + int * info); #endif // __APPLE__ namespace modmesh @@ -77,10 +129,53 @@ namespace modmesh #ifdef __APPLE__ using lapack_int_t = __LAPACK_int; +// Accelerate's new-LAPACK complex element types for CGEEV/ZGEEV. +using lapack_complex_float_t = __LAPACK_float_complex; +using lapack_complex_double_t = __LAPACK_double_complex; #else using lapack_int_t = int; +// On Linux the complex symbols above take void*; cast targets are void. +using lapack_complex_float_t = void; +using lapack_complex_double_t = void; #endif +namespace detail +{ + +/** + * Overloaded thin wrappers around the LAPACK *GEEV entry points, presenting a + * uniform call shape to EigenSystem. modmesh's Complex is + * layout-compatible with the LAPACK complex types, so the buffers are + * reinterpret_cast to the backend pointer type. + */ +inline void lapack_geev( + char jobvl, char jobvr, lapack_int_t n, float * a, lapack_int_t lda, float * wr, float * wi, float * vl, lapack_int_t ldvl, float * vr, lapack_int_t ldvr, float * work, lapack_int_t lwork, lapack_int_t * info) +{ + sgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, info); +} + +inline void lapack_geev( + char jobvl, char jobvr, lapack_int_t n, double * a, lapack_int_t lda, double * wr, double * wi, double * vl, lapack_int_t ldvl, double * vr, lapack_int_t ldvr, double * work, lapack_int_t lwork, lapack_int_t * info) +{ + dgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, info); +} + +inline void lapack_geev( + char jobvl, char jobvr, lapack_int_t n, Complex * a, lapack_int_t lda, Complex * w, Complex * vl, lapack_int_t ldvl, Complex * vr, lapack_int_t ldvr, Complex * work, lapack_int_t lwork, float * rwork, lapack_int_t * info) +{ + cgeev_( + &jobvl, &jobvr, &n, reinterpret_cast(a), &lda, reinterpret_cast(w), reinterpret_cast(vl), &ldvl, reinterpret_cast(vr), &ldvr, reinterpret_cast(work), &lwork, rwork, info); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast) +} + +inline void lapack_geev( + char jobvl, char jobvr, lapack_int_t n, Complex * a, lapack_int_t lda, Complex * w, Complex * vl, lapack_int_t ldvl, Complex * vr, lapack_int_t ldvr, Complex * work, lapack_int_t lwork, double * rwork, lapack_int_t * info) +{ + zgeev_( + &jobvl, &jobvr, &n, reinterpret_cast(a), &lda, reinterpret_cast(w), reinterpret_cast(vl), &ldvl, reinterpret_cast(vr), &ldvr, reinterpret_cast(work), &lwork, rwork, info); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast) +} + +} /* end namespace detail */ + } /* end namespace modmesh */ // vim: set ff=unix fenc=utf8 et sw=4 ts=4 sts=4: diff --git a/cpp/modmesh/linalg/pymod/wrap_EigenSystem.cpp b/cpp/modmesh/linalg/pymod/wrap_EigenSystem.cpp index 50b8fbdad..fe2fe5811 100644 --- a/cpp/modmesh/linalg/pymod/wrap_EigenSystem.cpp +++ b/cpp/modmesh/linalg/pymod/wrap_EigenSystem.cpp @@ -38,13 +38,14 @@ namespace python #ifdef MM_HAS_VENDOR_LAPACK +template class MODMESH_PYTHON_WRAPPER_VISIBILITY WrapEigenSystem - : public WrapBase + : public WrapBase, EigenSystem, std::shared_ptr>> { - using base_type = WrapBase; + using base_type = WrapBase, EigenSystem, std::shared_ptr>>; using wrapped_type = typename base_type::wrapped_type; - using array_type = SimpleArray; + using array_type = SimpleArray; friend base_type; @@ -52,7 +53,8 @@ class MODMESH_PYTHON_WRAPPER_VISIBILITY WrapEigenSystem }; /* end class WrapEigenSystem */ -WrapEigenSystem::WrapEigenSystem(pybind11::module & mod, char const * pyname, char const * pydoc) +template +WrapEigenSystem::WrapEigenSystem(pybind11::module & mod, char const * pyname, char const * pydoc) : base_type(mod, pyname, pydoc) { namespace py = pybind11; @@ -62,7 +64,7 @@ WrapEigenSystem::WrapEigenSystem(pybind11::module & mod, char const * pyname, ch py::init( [](array_type const & a, bool do_vl, bool do_vr) { - return std::make_unique(a, do_vl, do_vr); + return wrapped_type::construct(a, do_vl, do_vr); }), py::arg("a"), py::arg("do_vl") = true, @@ -70,24 +72,82 @@ WrapEigenSystem::WrapEigenSystem(pybind11::module & mod, char const * pyname, ch // Keep the input array's Python wrapper alive while the // EigenSystem lives, so m_matrix's C++ reference stays // valid for the lifetime of this instance. - py::keep_alive<1, 2>()) + py::keep_alive<1, 2>()); + + // Eigenvalues are exposed as real/imaginary parts for every element type; + // complex eigenvalues are split into wr/wi. + (*this) + .def_property_readonly("wr", &wrapped_type::wr) + .def_property_readonly("wi", &wrapped_type::wi); + + (*this) + .def("run", &wrapped_type::run) + .def_property_readonly( + "matrix", + &wrapped_type::matrix, + pybind11::return_value_policy::reference_internal) + .def_property_readonly( + "vl", + [](wrapped_type const & self) -> array_type const & + { + return self.vl(); + }) + .def_property_readonly( + "vr", + [](wrapped_type const & self) -> array_type const & + { + return self.vr(); + }) + .def( + "get_vl", + &wrapped_type::vl, + py::arg("suppress_exception") = false, + "Left eigenvectors; suppress_exception=True returns empty " + "instead of raising when do_vl=False.") + .def( + "get_vr", + &wrapped_type::vr, + py::arg("suppress_exception") = false, + "Right eigenvectors; suppress_exception=True returns empty " + "instead of raising when do_vr=False.") + .def_property_readonly("do_vl", &wrapped_type::do_vl) + .def_property_readonly("do_vr", &wrapped_type::do_vr) + .def_property_readonly("done", &wrapped_type::done); +} + +// Type-erased wrapper: constructs from a SimpleArrayPlex and dispatches on its +// runtime element type. Unlike WrapEigenSystem, it exposes wr/wi (real) and +// w (complex) together; the inapplicable ones raise at runtime. +class MODMESH_PYTHON_WRAPPER_VISIBILITY WrapEigenSystemPlex + : public WrapBase +{ + + using base_type = WrapBase; + using wrapped_type = typename base_type::wrapped_type; + + friend base_type; + + WrapEigenSystemPlex(pybind11::module & mod, char const * pyname, char const * pydoc); + +}; /* end class WrapEigenSystemPlex */ + +WrapEigenSystemPlex::WrapEigenSystemPlex(pybind11::module & mod, char const * pyname, char const * pydoc) + : base_type(mod, pyname, pydoc) +{ + namespace py = pybind11; + + (*this) .def( py::init( [](SimpleArrayPlex const & a, bool do_vl, bool do_vr) { - if (a.data_type() != DataType::Float64) - { - throw std::invalid_argument("EigenSystem: SimpleArray data type must be float64"); - } - // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast) - auto const * typed = reinterpret_cast(a.instance_ptr()); - return std::make_unique(*typed, do_vl, do_vr); + return std::make_unique(a, do_vl, do_vr); }), py::arg("a"), py::arg("do_vl") = true, py::arg("do_vr") = true, - // Keep the plex's Python wrapper alive while the EigenSystem - // lives; m_matrix references the array owned inside the plex. + // Keep the plex's Python wrapper alive while the EigenSystemPlex + // lives; the dispatched solver references the array it owns. py::keep_alive<1, 2>()); (*this) @@ -100,16 +160,12 @@ WrapEigenSystem::WrapEigenSystem(pybind11::module & mod, char const * pyname, ch .def_property_readonly("wi", &wrapped_type::wi) .def_property_readonly( "vl", - [](wrapped_type const & self) -> array_type const & - { - return self.vl(); - }) + [](wrapped_type const & self) + { return self.vl(); }) .def_property_readonly( "vr", - [](wrapped_type const & self) -> array_type const & - { - return self.vr(); - }) + [](wrapped_type const & self) + { return self.vr(); }) .def( "get_vl", &wrapped_type::vl, @@ -132,8 +188,18 @@ WrapEigenSystem::WrapEigenSystem(pybind11::module & mod, char const * pyname, ch void wrap_EigenSystem(pybind11::module & mod) { #ifdef MM_HAS_VENDOR_LAPACK - WrapEigenSystem::commit(mod, "EigenSystem", "Eigen problem solver"); + WrapEigenSystem::commit(mod, "EigenSystemFloat32", "Eigen problem solver (float32)"); + WrapEigenSystem::commit(mod, "EigenSystemFloat64", "Eigen problem solver (float64)"); + WrapEigenSystem>::commit(mod, "EigenSystemComplex64", "Eigen problem solver (complex64)"); + WrapEigenSystem>::commit(mod, "EigenSystemComplex128", "Eigen problem solver (complex128)"); + // EigenSystem is the type-erased, general entry point that infers the + // element type from a SimpleArrayPlex (C++ class: EigenSystemPlex). + WrapEigenSystemPlex::commit(mod, "EigenSystem", "Type-erased eigen problem solver"); #else // MM_HAS_VENDOR_LAPACK + mod.attr("EigenSystemFloat32") = pybind11::none(); + mod.attr("EigenSystemFloat64") = pybind11::none(); + mod.attr("EigenSystemComplex64") = pybind11::none(); + mod.attr("EigenSystemComplex128") = pybind11::none(); mod.attr("EigenSystem") = pybind11::none(); #endif // MM_HAS_VENDOR_LAPACK } diff --git a/modmesh/core.py b/modmesh/core.py index 55df05506..002ed253f 100644 --- a/modmesh/core.py +++ b/modmesh/core.py @@ -128,6 +128,10 @@ 'lu_solve', 'lu_inv', 'EigenSystem', + 'EigenSystemFloat32', + 'EigenSystemFloat64', + 'EigenSystemComplex64', + 'EigenSystemComplex128', 'KalmanStateInfoFp32', 'KalmanStateInfoFp64', 'KalmanStateInfoComplex64', diff --git a/tests/test_linalg_eigen.py b/tests/test_linalg_eigen.py index dc43f99b5..ab7363473 100644 --- a/tests/test_linalg_eigen.py +++ b/tests/test_linalg_eigen.py @@ -7,40 +7,198 @@ @unittest.skipIf(mm.EigenSystem is None, "mm.EigenSystem is not built (no vendor LAPACK)") -class TestEigenSystemTC(unittest.TestCase): - """Verify EigenSystem against DGEEV reference outputs. - - EigenSystem(A).run() populates wr/wi (real and imaginary parts of the - eigenvalues) and vl/vr (left/right eigenvector matrices) in DGEEV's - column-major layout (j-th column == j-th eigenvector). For a complex - conjugate pair (wi[j] > 0, wi[j+1] = -wi[j]) the j-th and (j+1)-th - columns hold the real and imaginary parts of the eigenvector; the - (j+1)-th eigenvector is the complex conjugate of the j-th. +class TestEigenSystemPlexTC(unittest.TestCase): + """Verify the type-erased EigenSystem surrogate (C++ EigenSystemPlex). + + mm.EigenSystem accepts a SimpleArrayPlex and dispatches on its runtime + element type to the matching typed EigenSystem. Eigenvalues are + exposed as real/imaginary parts via wr/wi for every element type. + """ + + REAL = ("float32", "float64") + COMPLEX = ("complex64", "complex128") + ALL = REAL + COMPLEX + + def _typed(self, arr, dtype): + # The matching typed (non-plex) solver, used as the oracle. + table = { + "float32": (mm.EigenSystemFloat32, mm.SimpleArrayFloat32), + "float64": (mm.EigenSystemFloat64, mm.SimpleArrayFloat64), + "complex64": (mm.EigenSystemComplex64, mm.SimpleArrayComplex64), + "complex128": (mm.EigenSystemComplex128, mm.SimpleArrayComplex128), + } + eig_cls, arr_cls = table[dtype] + return eig_cls(arr_cls(array=arr)) + + def test_dispatch_matches_typed_construction(self): + # For every supported dtype, the plex path must reproduce the typed + # EigenSystem path bit-for-bit: same dispatch, same GEEV call. + A_np = np.array([ + [2.0, 1.0, 0.5], + [0.0, 3.0, -1.0], + [0.0, 0.0, 5.0], + ]) + for dtype in self.ALL: + with self.subTest(dtype=dtype): + arr = np.ascontiguousarray(A_np, dtype=dtype) + solver = mm.EigenSystem(mm.SimpleArray(arr)) + self.assertFalse(solver.done) + solver.run() + self.assertTrue(solver.done) + typed = self._typed(arr, dtype) + typed.run() + np.testing.assert_array_equal( + np.array(solver.wr), np.array(typed.wr)) + np.testing.assert_array_equal( + np.array(solver.wi), np.array(typed.wi)) + np.testing.assert_array_equal( + np.array(solver.vl), np.array(typed.vl)) + np.testing.assert_array_equal( + np.array(solver.vr), np.array(typed.vr)) + + def test_rejects_unsupported_dtype(self): + # Integer element types have no GEEV; they must raise ValueError + # (std::invalid_argument) rather than dispatch. + for dtype in ("int32", "int64", "uint8"): + with self.subTest(dtype=dtype): + A = mm.SimpleArray(np.eye(2, dtype=dtype)) + with self.assertRaisesRegex( + ValueError, r"data type must be"): + mm.EigenSystem(A) + + def test_non_square_rejected(self): + # A supported dtype dispatches, then hits the square-2D guard inside + # the typed EigenSystem constructor. + A = mm.SimpleArray(np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])) + with self.assertRaisesRegex( + ValueError, r"must be a square 2D SimpleArray"): + mm.EigenSystem(A) + + def test_forwards_do_vl_do_vr_flags(self): + # do_vl/do_vr pass through to the dispatched solver unchanged. + A_np = np.array([[2.0, 0.0], [0.0, 3.0]]) + solver = mm.EigenSystem(mm.SimpleArray(A_np)) + self.assertTrue(solver.do_vl) + self.assertTrue(solver.do_vr) + + solver = mm.EigenSystem(mm.SimpleArray(A_np), do_vr=False) + solver.run() + self.assertTrue(solver.do_vl) + self.assertFalse(solver.do_vr) + with self.assertRaisesRegex( + RuntimeError, r"right eigenvectors were not computed"): + solver.vr # noqa: B018 + + def test_matrix_property_survives_input_gc(self): + # keep_alive<1, 2>() must keep the plex (and the typed array it owns) + # alive after the Python-side plex reference is dropped. + A_np = np.array([[3.0, 1.0], [0.0, 2.0]]) + plex = mm.SimpleArray(A_np) + solver = mm.EigenSystem(plex) + del plex + np.testing.assert_array_equal(np.array(solver.matrix), A_np) + + +class EigenSystemTB: + """Verify EigenSystem against *GEEV reference outputs. + + Subclasses bind a concrete element type by setting ``array_cls``, + ``eig_cls``, ``np_dtype``, ``is_complex`` and the comparison tolerances. + Real solvers (SGEEV/DGEEV) report eigenvalues as a real/imaginary split + (wr/wi) and pack a complex-conjugate eigenvector pair into two consecutive + real columns; complex solvers (CGEEV/ZGEEV) report a single complex w and + store eigenvectors directly as complex columns. The helpers below + normalize both layouts to complex arrays so each test body is + dtype-agnostic. """ - def _solve(self, A_np): - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A) + array_cls = None + eig_cls = None + np_dtype = None + is_complex = False + rtol = 1e-10 + atol = 1e-12 + + @classmethod + def get_complex_type(cls): + _ = {'float32': 'complex64', 'float64': 'complex128'} + name = np.dtype(cls.np_dtype).name + return _.get(name, name) + + def _array(self, A_np): + return self.array_cls( + array=np.ascontiguousarray(A_np, dtype=self.np_dtype)) + + def _solve(self, A_np, **kwargs): + solver = self.eig_cls(self._array(A_np), **kwargs) self.assertFalse(solver.done) solver.run() self.assertTrue(solver.done) - return (np.array(solver.wr), np.array(solver.wi), - np.array(solver.vl), np.array(solver.vr)) + return solver + + def _eigvals(self, solver): + # Eigenvalues as a complex ndarray: wr/wi for every element type. + wr = np.asarray(solver.wr, dtype=self.np_dtype) + wi = np.asarray(solver.wi, dtype=self.np_dtype) + return wr + 1j * wi + + def _columns_to_complex(self, mat, w): + # Normalize an eigenvector matrix to complex columns. For complex + # element types the columns are already complex; for real types *GEEV + # packs a conjugate pair (wi[j] > 0, wi[j+1] < 0) as + # v_j = M[:, j] + i M[:, j+1] and v_{j+1} = M[:, j] - i M[:, j+1]. + mat = np.asarray(mat) + if self.is_complex: + return mat.astype(self.np_dtype) + wi = np.imag(w) + n = mat.shape[0] + out = np.zeros((n, n), dtype=self.get_complex_type()) + j = 0 + while j < n: + if wi[j] == 0.0: + out[:, j] = mat[:, j] + j += 1 + else: + out[:, j] = mat[:, j] + 1j * mat[:, j + 1] + out[:, j + 1] = mat[:, j] - 1j * mat[:, j + 1] + j += 2 + return out + + def _assert_eigvals(self, got, expected): + # Compare eigenvalues as a multiset: order is unspecified and a sort + # key on the real/imag parts is fragile near ties (e.g. +/- i, where + # rounding noise in the zero real part flips the order). Greedily + # match each expected value to its nearest computed one instead. + got = list(np.asarray(got, dtype=self.get_complex_type()).ravel()) + expected = np.asarray(expected, dtype=self.get_complex_type()).ravel() + self.assertEqual(len(got), expected.size) + for e in expected: + diffs = np.abs(np.asarray(got) - e) + k = int(np.argmin(diffs)) + tol = self.atol + self.rtol * abs(e) + self.assertLessEqual( + diffs[k], tol, + msg=f"no computed eigenvalue near {e} (closest off by " + f"{diffs[k]:.3e})") + got.pop(k) def test_diagonal_eigenvalues_and_basis_eigenvectors(self): # diag(1, 2, 3, 4): eigenvalues are the diagonal; eigenvectors are - # standard basis vectors up to permutation/sign. + # standard basis vectors up to permutation/phase. A_np = np.diag([1.0, 2.0, 3.0, 4.0]) - wr, wi, _vl, vr = self._solve(A_np) - np.testing.assert_allclose(np.sort(wr), [1.0, 2.0, 3.0, 4.0], - rtol=1e-12, atol=1e-14) - np.testing.assert_allclose(wi, np.zeros(4), atol=1e-14) - # Each column must be a (signed) standard basis vector. + solver = self._solve(A_np) + w = self._eigvals(solver) + self._assert_eigvals(w, [1.0, 2.0, 3.0, 4.0]) + np.testing.assert_allclose(np.imag(w), np.zeros(4), + atol=10 * self.atol) + vr = self._columns_to_complex(solver.vr, w) + # Each column is a (phased) standard basis vector: one unit-magnitude + # entry and the rest zero. for j in range(4): col = vr[:, j] np.testing.assert_allclose(np.abs(col).sum(), 1.0, - rtol=1e-12, atol=1e-14) - self.assertEqual(np.count_nonzero(np.abs(col) > 1e-12), 1) + rtol=self.rtol, atol=10 * self.atol) + self.assertEqual(np.count_nonzero(np.abs(col) > 1e-4), 1) def test_symmetric_known_spectrum_reconstructs_eigenpairs(self): # The fixture Q below is frozen (not generated at runtime) so the test @@ -71,119 +229,94 @@ def test_symmetric_known_spectrum_reconstructs_eigenpairs(self): # Check the fixture is orthogonal first. np.testing.assert_allclose(Q @ Q.T, np.eye(n), atol=1e-14) A_np = Q @ np.diag(lam) @ Q.T - wr, wi, _vl, vr = self._solve(A_np) - np.testing.assert_allclose(np.sort(wr), np.sort(lam), - rtol=1e-10, atol=1e-12) - np.testing.assert_allclose(wi, np.zeros(n), atol=1e-12) + solver = self._solve(A_np) + w = self._eigvals(solver) + self._assert_eigvals(w, lam) + np.testing.assert_allclose(np.imag(w), np.zeros(n), + atol=10 * self.atol) # Per-column eigenpair: A v_j = lambda_j v_j. + Ac = A_np.astype(complex) + vr = self._columns_to_complex(solver.vr, w) for j in range(n): - np.testing.assert_allclose(A_np @ vr[:, j], wr[j] * vr[:, j], - rtol=1e-10, atol=1e-12) + np.testing.assert_allclose(Ac @ vr[:, j], w[j] * vr[:, j], + rtol=self.rtol, atol=10 * self.atol) def test_2x2_rotation_complex_conjugate_pair(self): - # 2x2 90-degree rotation matrix has eigenvalues +/- i. This fixture - # exercises DGEEV's packing of complex conjugate eigenvectors into - # consecutive real columns. + # 2x2 90-degree rotation matrix has eigenvalues +/- i. For real + # element types this exercises *GEEV's packing of a complex conjugate + # eigenvector pair into consecutive real columns; for complex types the + # eigenvectors come back complex directly. A_np = np.array([[0.0, -1.0], [1.0, 0.0]], dtype='float64') - wr, wi, _vl, vr = self._solve(A_np) - np.testing.assert_allclose(wr, np.zeros(2, dtype='float64'), - atol=1e-14) - np.testing.assert_allclose(np.sort(wi), [-1.0, 1.0], atol=1e-14) - j = int(np.argmax(wi)) # column with positive imaginary part - self.assertGreater(wi[j], 0.0) - self.assertAlmostEqual(wi[j] + wi[j + 1], 0.0, places=14) - v = vr[:, j] + 1j * vr[:, j + 1] - lam = wr[j] + 1j * wi[j] - np.testing.assert_allclose(A_np @ v, lam * v, rtol=1e-12, atol=1e-14) - # The (j+1)-th eigenvector is the conjugate of the j-th and corresponds - # to the conjugate eigenvalue. - v_conj = vr[:, j] - 1j * vr[:, j + 1] - lam_conj = wr[j] - 1j * wi[j] - np.testing.assert_allclose(A_np @ v_conj, lam_conj * v_conj, - rtol=1e-12, atol=1e-14) + solver = self._solve(A_np) + w = self._eigvals(solver) + self._assert_eigvals(w, [1j, -1j]) + Ac = A_np.astype(complex) + vr = self._columns_to_complex(solver.vr, w) + for j in range(2): + np.testing.assert_allclose(Ac @ vr[:, j], w[j] * vr[:, j], + rtol=self.rtol, atol=10 * self.atol) def test_left_eigenvectors_satisfy_left_equation(self): - # Left eigenvectors u_j (column vector) satisfy u_j^T A = lambda_j - # u_j^T (i.e., A^T u_j = lambda_j u_j) for real eigenvalues (general - # case is A^H u_j = conj(lambda_j) u_j). Use a non-symmetric matrix - # with all-real eigenvalues: upper-triangular -> spectrum is the - # diagonal. + # Left eigenvectors u_j satisfy A^H u_j = conj(lambda_j) u_j. Use a + # non-symmetric upper-triangular matrix whose spectrum is its diagonal. A_np = np.array([ [2.0, 1.0, 0.5], [0.0, 3.0, -1.0], [0.0, 0.0, 5.0], ], dtype='float64') - wr, wi, vl, _vr = self._solve(A_np) - np.testing.assert_allclose(np.sort(wr), [2.0, 3.0, 5.0], - rtol=1e-12, atol=1e-12) - np.testing.assert_allclose(wi, np.zeros(3, dtype='float64'), - atol=1e-12) + solver = self._solve(A_np) + w = self._eigvals(solver) + self._assert_eigvals(w, [2.0, 3.0, 5.0]) + Ah = A_np.astype(complex).conj().T + vl = self._columns_to_complex(solver.vl, w) for j in range(3): - np.testing.assert_allclose(A_np.T @ vl[:, j], - wr[j] * vl[:, j], - rtol=1e-10, atol=1e-12) + np.testing.assert_allclose(Ah @ vl[:, j], + np.conj(w[j]) * vl[:, j], + rtol=self.rtol, atol=10 * self.atol) def test_rejects_non_square_and_non_2d_inputs(self): - # Non-square 2D, 1D, and 3D inputs must all raise the same shape - # error from the EigenSystem constructor. - A_rect = mm.SimpleArrayFloat64(array=np.array( - [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype="float64")) - with self.assertRaisesRegex( - ValueError, r"must be a square 2D SimpleArray"): - mm.EigenSystem(A_rect) - - A_1d = mm.SimpleArrayFloat64(array=np.array( - [1.0, 2.0, 3.0], dtype="float64")) - with self.assertRaisesRegex( - ValueError, r"must be a square 2D SimpleArray"): - mm.EigenSystem(A_1d) - - A_3d = mm.SimpleArrayFloat64(array=np.zeros((2, 2, 2), - dtype="float64")) - with self.assertRaisesRegex( - ValueError, r"must be a square 2D SimpleArray"): - mm.EigenSystem(A_3d) + # Non-square 2D, 1D, and 3D inputs must all raise the same shape error + # from the EigenSystem constructor. + for bad in (np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]), + np.array([1.0, 2.0, 3.0]), + np.zeros((2, 2, 2))): + A = self._array(bad) + with self.assertRaisesRegex( + ValueError, r"must be a square 2D SimpleArray"): + self.eig_cls(A) def test_accessors_before_run_are_inert(self): - # Construct but do not call run(); done must be False and wr must be - # finite-valued (DGEEV hasn't written), documenting that run() is + # Construct but do not call run(); done must be False and the + # eigenvalue accessor must be finite, documenting that run() is # required before reading results. - A_np = np.diag([1.0, 2.0]) - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A) + A = self._array(np.diag([1.0, 2.0])) + solver = self.eig_cls(A) self.assertFalse(solver.done) - wr = solver.wr.ndarray - self.assertTrue(np.all(np.isfinite(wr))) + self.assertTrue(np.all(np.isfinite(self._eigvals(solver)))) def test_matrix_property_survives_input_gc(self): - # Construct solver, then drop the Python-side reference to A. - # solver.matrix must still equal the original array, confirming that - # keep_alive<1, 2>() on the constructor prevents the C++ m_matrix - # reference from dangling. - A_np = np.array([[3.0, 1.0], [0.0, 2.0]]) - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A) + # Drop the Python-side reference to A; solver.matrix must still equal + # the original, confirming keep_alive<1, 2>() keeps m_matrix valid. + A_np = np.ascontiguousarray(np.array([[3.0, 1.0], [0.0, 2.0]]), + dtype=self.np_dtype) + A = self.array_cls(array=A_np) + solver = self.eig_cls(A) del A np.testing.assert_array_equal(solver.matrix.ndarray, A_np) def test_default_constructor_flags_both_true(self): # Guard against an accidental default-flag flip. - A = mm.SimpleArrayFloat64(array=np.diag([1.0, 2.0])) - solver = mm.EigenSystem(A) + solver = self.eig_cls(self._array(np.diag([1.0, 2.0]))) self.assertTrue(solver.do_vl) self.assertTrue(solver.do_vr) def test_skip_vr_does_not_compute_right_eigenvectors(self): - # do_vr=False must keep wr/wi/vl intact and reject both vr accessors. - A_np = np.array([[2.0, 0.0], [0.0, 3.0]], dtype="float64") - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A, do_vr=False) - solver.run() - self.assertTrue(solver.done) + # do_vr=False must keep eigenvalues/vl intact and reject both vr + # accessors. + solver = self._solve(np.array([[2.0, 0.0], [0.0, 3.0]]), do_vr=False) self.assertTrue(solver.do_vl) self.assertFalse(solver.do_vr) - np.testing.assert_allclose(np.sort(np.array(solver.wr)), - [2.0, 3.0], atol=1e-14) + self._assert_eigvals(self._eigvals(solver), [2.0, 3.0]) self.assertEqual(np.array(solver.vl).shape, (2, 2)) self.assertEqual(np.array(solver.get_vl()).shape, (2, 2)) with self.assertRaisesRegex( @@ -195,18 +328,16 @@ def test_skip_vr_does_not_compute_right_eigenvectors(self): def test_skip_vl_does_not_compute_left_eigenvectors(self): # Mirror of test_skip_vr_*; checks A v = lambda v stays exact. - A_np = np.array([[2.0, 1.0], [0.0, 3.0]], dtype="float64") - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A, do_vl=False) - solver.run() - self.assertTrue(solver.done) + A_np = np.array([[2.0, 1.0], [0.0, 3.0]], dtype='float64') + solver = self._solve(A_np, do_vl=False) self.assertFalse(solver.do_vl) self.assertTrue(solver.do_vr) - wr = np.array(solver.wr) - vr = np.array(solver.vr) + w = self._eigvals(solver) + Ac = A_np.astype(complex) + vr = self._columns_to_complex(solver.vr, w) for j in range(2): - np.testing.assert_allclose(A_np @ vr[:, j], wr[j] * vr[:, j], - rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(Ac @ vr[:, j], w[j] * vr[:, j], + rtol=self.rtol, atol=10 * self.atol) with self.assertRaisesRegex( RuntimeError, r"left eigenvectors were not computed"): solver.vl # noqa: B018 @@ -215,13 +346,10 @@ def test_skip_vl_does_not_compute_left_eigenvectors(self): solver.get_vl() def test_skip_both_computes_only_eigenvalues(self): - # Exercises the 3*n workspace path (both jobvl=jobvr='N'). - A_np = np.diag([1.0, 5.0, 9.0]) - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A, do_vl=False, do_vr=False) - solver.run() - np.testing.assert_allclose(np.sort(np.array(solver.wr)), - [1.0, 5.0, 9.0], atol=1e-14) + # Exercises the eigenvalue-only workspace path (jobvl=jobvr='N'). + solver = self._solve(np.diag([1.0, 5.0, 9.0]), + do_vl=False, do_vr=False) + self._assert_eigvals(self._eigvals(solver), [1.0, 5.0, 9.0]) with self.assertRaises(RuntimeError): solver.vl # noqa: B018 with self.assertRaises(RuntimeError): @@ -232,12 +360,9 @@ def test_skip_both_computes_only_eigenvalues(self): solver.get_vr() def test_get_methods_match_property_when_computed(self): - # Property and get_v* must alias the same matrix; suppress flag - # is a no-op when the matrix exists. - A_np = np.diag([2.0, 4.0]) - A = mm.SimpleArrayFloat64(array=A_np) - solver = mm.EigenSystem(A) - solver.run() + # Property and get_v* must alias the same matrix; suppress flag is a + # no-op when the matrix exists. + solver = self._solve(np.diag([2.0, 4.0])) np.testing.assert_array_equal(np.array(solver.vl), np.array(solver.get_vl())) np.testing.assert_array_equal(np.array(solver.vr), @@ -251,9 +376,7 @@ def test_get_methods_match_property_when_computed(self): def test_get_methods_default_argument_raises(self): # Only suppress_exception=True silences the exception. - A = mm.SimpleArrayFloat64(array=np.diag([1.0, 2.0])) - solver = mm.EigenSystem(A, do_vl=False, do_vr=False) - solver.run() + solver = self._solve(np.diag([1.0, 2.0]), do_vl=False, do_vr=False) with self.assertRaises(RuntimeError): solver.get_vl() with self.assertRaises(RuntimeError): @@ -265,91 +388,54 @@ def test_get_methods_default_argument_raises(self): def test_suppress_exception_returns_empty_when_not_computed(self): # suppress_exception=True returns an empty placeholder, not data. - A = mm.SimpleArrayFloat64(array=np.diag([1.0, 2.0])) - solver = mm.EigenSystem(A, do_vl=False, do_vr=False) - solver.run() + solver = self._solve(np.diag([1.0, 2.0]), do_vl=False, do_vr=False) empty_vl = np.array(solver.get_vl(suppress_exception=True)) empty_vr = np.array(solver.get_vr(suppress_exception=True)) self.assertEqual(empty_vl.size, 0) self.assertEqual(empty_vr.size, 0) -@unittest.skipIf(mm.EigenSystem is None, - "mm.EigenSystem is not built (no vendor LAPACK)") -class TestEigenSystemPlexTC(unittest.TestCase): - """Verify EigenSystem construction from a type-erased SimpleArray - - The type-erased SimpleArray is SimpleArrayPlex in C++. - """ - - def test_plex_float64_matches_typed_construction(self): - # A float64 plex must construct EigenSystem and yield results identical - # to the typed SimpleArrayFloat64 path: same data, same DGEEV call. - A_np = np.array([ - [2.0, 1.0, 0.5], - [0.0, 3.0, -1.0], - [0.0, 0.0, 5.0], - ], dtype="float64") - plex = mm.SimpleArray(A_np) - solver = mm.EigenSystem(plex) - self.assertFalse(solver.done) - solver.run() - self.assertTrue(solver.done) - # Identical to the typed SimpleArrayFloat64 construction. - typed = mm.EigenSystem(mm.SimpleArrayFloat64(array=A_np)) - typed.run() - for name in ("wr", "wi", "vl", "vr"): - np.testing.assert_array_equal( - np.array(getattr(solver, name)), - np.array(getattr(typed, name))) - # Sanity: the eigenvalues are the diagonal entries. - np.testing.assert_allclose(np.sort(np.array(solver.wr)), - [2.0, 3.0, 5.0], rtol=1e-12, atol=1e-12) - - def test_plex_rejects_non_float64_dtype(self): - # The plex overload only accepts float64. float32 and integer element - # types must raise ValueError (std::invalid_argument). - for dtype in ("float32", "int32"): - A_np = np.array([[2.0, 0.0], [0.0, 3.0]], dtype=dtype) - A = mm.SimpleArray(A_np) - with self.assertRaisesRegex( - ValueError, r"data type must be float64"): - mm.EigenSystem(A) - - def test_plex_non_square_rejected(self): - # The float64 dtype check passes, so a non-square plex must still hit - # the square-2D guard in the EigenSystem constructor. - A_np = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype="float64") - A = mm.SimpleArray(A_np) - with self.assertRaisesRegex( - ValueError, r"must be a square 2D SimpleArray"): - mm.EigenSystem(A) - - def test_plex_forwards_do_vl_do_vr_flags(self): - # do_vl/do_vr must pass through the plex overload unchanged: the - # default keeps both true, and do_vr=False suppresses vr. - A_np = np.array([[2.0, 0.0], [0.0, 3.0]], dtype="float64") - solver = mm.EigenSystem(mm.SimpleArray(A_np)) - self.assertTrue(solver.do_vl) - self.assertTrue(solver.do_vr) - - solver = mm.EigenSystem(mm.SimpleArray(A_np), do_vr=False) - solver.run() - self.assertTrue(solver.do_vl) - self.assertFalse(solver.do_vr) - with self.assertRaisesRegex( - RuntimeError, r"right eigenvectors were not computed"): - solver.vr # noqa: B018 - - def test_plex_matrix_property_survives_input_gc(self): - # Mirror of test_matrix_property_survives_input_gc for the plex - # overload: m_matrix references the array owned inside the plex, so - # py::keep_alive<1, 2>() must stop it from dangling once the - # Python-side plex reference is dropped. - A_np = np.array([[3.0, 1.0], [0.0, 2.0]], dtype="float64") - plex = mm.SimpleArray(A_np) - solver = mm.EigenSystem(plex) - del plex - np.testing.assert_array_equal(solver.matrix.ndarray, A_np) +@unittest.skipIf(mm.EigenSystemFloat32 is None, + "EigenSystem is not built (no vendor LAPACK)") +class TestLinalgEigenSystemFloat32TC(EigenSystemTB, unittest.TestCase): + array_cls = mm.SimpleArrayFloat32 + eig_cls = mm.EigenSystemFloat32 + np_dtype = np.float32 + is_complex = False + rtol = 1e-4 + atol = 1e-5 + + +@unittest.skipIf(mm.EigenSystemFloat64 is None, + "EigenSystem is not built (no vendor LAPACK)") +class TestLinalgEigenSystemFloat64TC(EigenSystemTB, unittest.TestCase): + array_cls = mm.SimpleArrayFloat64 + eig_cls = mm.EigenSystemFloat64 + np_dtype = np.float64 + is_complex = False + rtol = 1e-10 + atol = 1e-12 + + +@unittest.skipIf(mm.EigenSystemComplex64 is None, + "EigenSystem is not built (no vendor LAPACK)") +class TestLinalgEigenSystemComplex64TC(EigenSystemTB, unittest.TestCase): + array_cls = mm.SimpleArrayComplex64 + eig_cls = mm.EigenSystemComplex64 + np_dtype = np.complex64 + is_complex = True + rtol = 1e-4 + atol = 1e-5 + + +@unittest.skipIf(mm.EigenSystemComplex128 is None, + "EigenSystem is not built (no vendor LAPACK)") +class TestLinalgEigenSystemComplex128TC(EigenSystemTB, unittest.TestCase): + array_cls = mm.SimpleArrayComplex128 + eig_cls = mm.EigenSystemComplex128 + np_dtype = np.complex128 + is_complex = True + rtol = 1e-10 + atol = 1e-12 # vim: set ff=unix fenc=utf8 et sw=4 ts=4 sts=4: