Skip to content

Add EigenSystem with multi-dtype support and shared_ptr management#833

Merged
yungyuc merged 1 commit into
solvcon:masterfrom
yungyuc:feat/eigentmpl
May 26, 2026
Merged

Add EigenSystem with multi-dtype support and shared_ptr management#833
yungyuc merged 1 commit into
solvcon:masterfrom
yungyuc:feat/eigentmpl

Conversation

@yungyuc
Copy link
Copy Markdown
Member

@yungyuc yungyuc commented May 26, 2026

Add float32/complex64/complex128 variants of EigenSystem alongside the existing float64 path. Complex eigenvalues are returned as separate wr/wi real arrays, and the Plex wrapper is renamed to EigenSystem in Python. The class is now managed exclusively via std::shared_ptr through a construct() factory function,

The type-erased EigenSystemPlex implementation is consolidated into EigenSystem.cpp to save compilation time.

Add more LAPACK type compatbility in lapack_compat.hpp and update the Python tests accordingly.

For issue #790.

Add float32/complex64/complex128 variants of EigenSystem alongside the existing
float64 path.  Complex eigenvalues are returned as separate wr/wi real arrays,
and the Plex wrapper is renamed to `EigenSystem` in Python.  The class is now
managed exclusively via std::shared_ptr through a construct() factory function,

The type-erased EigenSystemPlex implementation is consolidated into
`EigenSystem.cpp` to save compilation time.

Add more LAPACK type compatbility in `lapack_compat.hpp` and update the Python
tests accordingly.
@yungyuc yungyuc self-assigned this May 26, 2026
@yungyuc yungyuc added the array Multi-dimensional array implementation label May 26, 2026
@yungyuc yungyuc moved this from Todo to In Progress in tensor operations May 26, 2026
// Reinterpret the plex's type-erased instance as the typed SimpleArray<T> it
// actually holds; the caller must have verified data_type() == T.
template <typename T>
static SimpleArray<T> const & typed_array(SimpleArrayPlex const & plex)
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Static function template to support type-erased SimpleArray.

static constexpr bool is_complex = is_complex_v<T>;

explicit EigenSystem(array_type const & matrix, bool do_vl = true, bool do_vr = true);
static std::shared_ptr<EigenSystem> construct(array_type const & matrix, bool do_vl = true, bool do_vr = true)
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Factory function template to make EigenSystem exclusively managed by using std::shared_ptr.

throw std::runtime_error(std::format(
"EigenSystem::run: DGEEV workspace query failed with info={}",
static_cast<int64_t>(info)));
// CGEEV/ZGEEV return eigenvalues as a single complex array and need a
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[CZ]GEEV require a single complex array as the working space. To keep the external interface the same, the implementation here differs between complex and scalar input.


/**
* Overloaded thin wrappers around the LAPACK *GEEV entry points, presenting a
* uniform call shape to EigenSystem<T>. modmesh's Complex<T> is
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Complex<T> uses the same memory layout as std::complex.

WrapEigenSystem<Complex<double>>::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");
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Type-erased EigenSystemPlex is wrapped to Python as EigenSystem.

@yungyuc yungyuc marked this pull request as ready for review May 26, 2026 13:25
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):
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Test the type-erased EigenSystem.


REAL = ("float32", "float64")
COMPLEX = ("complex64", "complex128")
ALL = REAL + COMPLEX
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All 4 types of arrays.


@classmethod
def get_complex_type(cls):
_ = {'float32': 'complex64', 'float64': 'complex128'}
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use corresponding complex types.

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):
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Test for all 4 types with the explicitly typed EigenSystem* classes.

@yungyuc yungyuc changed the title Add EigenSystem with multi-dtype support and shared_ptr management Add EigenSystem with multi-dtype support and shared_ptr management May 26, 2026
@yungyuc yungyuc merged commit a0cbe6a into solvcon:master May 26, 2026
41 of 49 checks passed
@github-project-automation github-project-automation Bot moved this from In Progress to Done in tensor operations May 26, 2026
@yungyuc yungyuc deleted the feat/eigentmpl branch May 26, 2026 14:02
@yungyuc yungyuc linked an issue May 26, 2026 that may be closed by this pull request
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

array Multi-dimensional array implementation

Projects

Status: Done

Development

Successfully merging this pull request may close these issues.

Add float32, complex64, complex128 variants of EigenSystem

1 participant