Add LU decomposition with partial pivoting by tigercosmos · Pull Request #735 · solvcon/modmesh

tigercosmos · 2026-04-18T18:28:48Z

Summary

Implement LU factorization with partial pivoting (PA=LU) for general (non-symmetric) matrices in cpp/modmesh/linalg/lu_factorization.hpp
Add lu_factorization(), lu_solve(), and lu_inv() free functions with Python bindings for float32, float64, complex64, and complex128 types
Add .solve(b) and .inv() convenience methods on SimpleArray for floating-point and complex types (integer arrays excluded)

Test plan

test_lu_solve_and_inv_dtypes: parametric test across all 4 supported dtypes (float64, float32, complex128, complex64) verifying solve and inverse
test_lu_solve_2d_rhs: multi-RHS solve with 5×5 matrix and 5×4 RHS (n,m both > 3)
test_lu_singular: singular matrix raises RuntimeError
test_lu_invalid_inputs: non-square, 1D, dimension mismatch, wrong rank all raise ValueError
test_simplearray_int_no_solve_inv: integer SimpleArray does not expose solve/inv
All 35 tests pass (pytest tests/test_linalg.py)

Document

Please also check the Markdown file (gist) for detailed explanation of the math and code.

The PR is part of #719.

tigercosmos · 2026-04-18T18:49:19Z

+/**
+ * LU decomposition with partial pivoting for general (non-symmetric) matrices.
+ *
+ * Given an n-by-n matrix A, this class computes the factorization PA = LU,
+ * where:
+ *   - P is a permutation matrix (represented as a pivot index vector),
+ *   - L is a unit lower triangular matrix (diagonal elements are implicitly 1),
+ *   - U is an upper triangular matrix.
+ *
+ * L and U are stored compactly in a single n-by-n array: the strictly lower
+ * triangle holds L (excluding the unit diagonal), and the upper triangle
+ * (including the diagonal) holds U.
+ *
+ * The pivot vector piv[k] records that row k was swapped with row piv[k]
+ * during the k-th elimination step.  Swaps are applied sequentially from
+ * k = 0 to k = n-1.
+ *
+ * Supported element types: float, double, Complex<float>, Complex<double>.
+ */


I tried to make the comment clear. The PR description also includes a HTML explanation file. Almost all code are generated by AI, but I have checked all logic.

tigercosmos · 2026-04-18T18:50:45Z

+     * @throws std::invalid_argument if a is not a square 2D array.
+     * @throws std::runtime_error    if the matrix is singular or near-singular.
+     */
+    static std::pair<array_type, std::vector<ssize_t>> factorize(array_type const & a);


tigercosmos · 2026-04-18T18:51:11Z

+     *         dimensions are incompatible.
+     * @throws std::runtime_error    if a is singular.
+     */
+    static array_type solve(array_type const & a, array_type const & b);


Solver's entry.

tigercosmos · 2026-04-18T18:53:01Z

+template <typename T>
+SimpleArray<T> lu_solve(SimpleArray<T> const & a, SimpleArray<T> const & b)
+{
+    return detail::Lu<T>::solve(a, b);
+}


Export the free function.

tigercosmos · 2026-04-18T18:53:39Z

+        cases = [
+            # (np_dtype, SimpleArray class, matrix size, rtol, atol)
+            (np.float64, mm.SimpleArrayFloat64, 4, 1e-10, 1e-10),
+            (np.float32, mm.SimpleArrayFloat32, 3, 1e-5, 1e-5),
+            (np.complex128, mm.SimpleArrayComplex128, 3, 1e-10, 1e-10),
+            (np.complex64, mm.SimpleArrayComplex64, 2, 1e-5, 1e-5),
+        ]


Try to cover different dtypes and matrix sizes.

tigercosmos · 2026-04-18T18:54:47Z

        np.testing.assert_allclose(ps_upd, ps_upd_np, atol=1e-12, rtol=0.0)

+
+class TestLuSolver(unittest.TestCase):


There could be more test cases. I intended keep the most important ones to keep the PR small.

In this case, just add more tests without worrying about the diff length. The tests are necessary.

Sure, added more test cases.

tigercosmos · 2026-04-18T19:02:11Z

@yungyuc Please take a look. The PR slightly exceeds 500 LOC, but some of the code are just document comments.

yungyuc

Points to address:

It's OK to add more testing code even though there are many lines already. The testing code is the key to the work.
Use explicit test fixture.
Clarify why testing for the include file existence.

yungyuc · 2026-04-19T01:08:12Z

        np.testing.assert_allclose(ps_upd, ps_upd_np, atol=1e-12, rtol=0.0)

+
+class TestLuSolver(unittest.TestCase):


In this case, just add more tests without worrying about the diff length. The tests are necessary.

yungyuc · 2026-04-19T01:09:45Z

 #include <modmesh/buffer/pymod/buffer_pymod.hpp> // Must be the first include.

 #include <modmesh/buffer/pymod/array_common.hpp>
+#if __has_include(<modmesh/linalg/lu_factorization.hpp>)


Why do you need to test for the existence of the include file? In what scenario it would not exist?

Indeed. Remove it.

I guess you meant that it is not necessary.

Yes, thank you.

yungyuc · 2026-04-19T02:24:51Z

Document

Please also check the HTML file (gist) for detailed explanation of the math and code.

It's not easy to read the raw HTML. Gist does not seem to allow rendering it.

Using markdown with mathjax (e.g., $$\mathrm{A} x = b$$) would be easier to read, but it may take additional time to write. We can discuss what is a better approach.

tigercosmos · 2026-04-19T05:25:15Z

Using markdown with mathjax would be easier to read, but it may take additional time to write. We can discuss what is a better approach.

I think gist will only show the raw file anyway? Maybe I am wrong.
I choose HTML files because it can have the most flexibility on visualization, such as charts, colors or even animations.

In this case, I would suggest you download the html file, and use a browser to open it.

yungyuc · 2026-04-19T05:34:03Z

Gist render mathjax in markdown. See https://gist.github.com/yungyuc/3f5fd600c91a536d881c581a2045eb06 for an example.

tigercosmos · 2026-04-19T05:47:41Z

@yungyuc Thanks, I will try using markdown next time.

yungyuc · 2026-04-19T05:51:12Z

No problem. Could you rewrite it in markdown this time? It is very inconvenient to download and review HTML files. It is also insecure.

tigercosmos · 2026-04-19T08:10:11Z

No problem. Could you rewrite it in markdown this time? It is very inconvenient to download and review HTML files. It is also insecure.

Sure, I have update the gist to Markdown.

tigercosmos · 2026-04-19T12:25:34Z

@yungyuc I have added more test cases in fixtures. Please take a look, thanks!

yungyuc

A point to address in the engine:

Lu::format_shape() is long enough to be put outside the class declaration of Lu.

Points to address for testing:

Always add dtype when creating ndarray.
Clarify what does "complex-path" mean.
Make sure the factorization error out correctly when the system has a very small eigenvalue (nearly singular).

yungyuc · 2026-04-20T13:59:42Z

+            [2.0, 1.0, 1.0],
+            [4.0, -6.0, 0.0],
+            [-2.0, 7.0, 2.0],
+        ])


Always add dtype='string' when creating ndarray. The explicit dtype specification makes maintenance and grepping easier.

Fixed, add dtype for all arrays.

yungyuc · 2026-04-20T14:02:06Z

+            [1.0, 1.0, 1.0],
+            [0.0, 1.0, 2.0],
+        ])
+        # 3x3 complex matrix for complex-path coverage.


What does "complex-path" mean? Please provide a reference in the PR comment.

Modified the comment. It should meant testing Lu<Complex<double>>.

yungyuc · 2026-04-20T14:06:06Z

+                ValueError, r"must be a square 2D SimpleArray"):
+            mm.lu_factorization(A_1d)
+
+    def test_factorize_rejects_singular_duplicate_row(self):


It's important to test for singular system. Duplicated rows or columns make the system singular, but very small eigenvalues also make the system almost singular.

Please add a test for tiny eigenvalue.

Added a new test_factorize_rejects_near_singular_tiny_eigenvalue.

yungyuc · 2026-04-20T14:07:53Z

 #include <modmesh/buffer/pymod/buffer_pymod.hpp> // Must be the first include.

 #include <modmesh/buffer/pymod/array_common.hpp>
+#if __has_include(<modmesh/linalg/lu_factorization.hpp>)


I guess you meant that it is not necessary.

yungyuc · 2026-04-20T14:09:39Z

+    using array_type = SimpleArray<value_type>;
+
+    // Helper: format a SimpleArray shape as "(d0, d1, ...)" for error messages.
+    static std::string format_shape(array_type const & arr)


Lu::format_shape() is long enough to be put outside the class declaration of Lu. But if it is intentional to put it here, it's short enough to be OK.

yungyuc · 2026-04-20T14:11:06Z

+{
+
+    static_assert(
+        std::is_floating_point_v<T> || is_complex_v<T>,


In the future we can consider to implement is_real_v by making it an alias to std::is_floating_point_v. Then the error message can simply write "... a real or complex number type".

tigercosmos

@yungyuc Please take a look. Thanks!

tigercosmos · 2026-04-22T18:22:30Z

 constexpr bool is_complex_v = is_complex<T>::value;

+template <typename T>
+constexpr bool is_real_v = is_real<T>::value;


Added is_real_v.

tigercosmos · 2026-04-22T18:23:04Z

+            [1.0, 1.0, 1.0],
+            [0.0, 1.0, 2.0],
+        ], dtype=np.float64)
+        # 3x3 complex matrix to exercise Lu<T> instantiated for


Fixed the comment.

tigercosmos · 2026-04-22T18:24:31Z

+                RuntimeError, r"singular or near-singular"):
+            mm.lu_factorization(A)
+
+    def test_factorize_rejects_near_singular_tiny_eigenvalue(self):


New added test case for the small eigenvalue.

Good addition.

tigercosmos · 2026-04-22T18:25:05Z

+            [5.0, 6.0, 8.0, 8.0],
+            [9.0, 10.0, 11.0, 16.0],
+            [13.0, 15.0, 16.0, 17.0],
+        ], dtype=np.float64)


Always added dtype.

It will help reader to always know what data type the array uses.

Please use the string version dtype='float64' (double quotes is OK if you like: dtype="float64") for fundamental types.

In terms of functionality it does not matter to use string or object for dtype. It's just for style consistency.

tigercosmos · 2026-04-22T18:32:28Z

+        // Absolute threshold (~100 * machine eps); works for well-scaled inputs.
+        // TODO: make it relative to matrix/column magnitude for better robustness.
+        real_type const singular_tol = real_type(100) * eps;
+        if (pivot <= singular_tol)


One more changed compared with the previous review.

yungyuc

Only a minor formatting issue remains. If you want to fix it in the future, I can merge now:

Use the string version for dtype of fundamental types, e.g., dtype='float64' (double quotes is OK if you like: dtype="float64") for fundamental types.

yungyuc · 2026-04-26T02:19:02Z

+                RuntimeError, r"singular or near-singular"):
+            mm.lu_factorization(A)
+
+    def test_factorize_rejects_near_singular_tiny_eigenvalue(self):


Good addition.

yungyuc · 2026-04-26T02:22:14Z

+            [5.0, 6.0, 8.0, 8.0],
+            [9.0, 10.0, 11.0, 16.0],
+            [13.0, 15.0, 16.0, 17.0],
+        ], dtype=np.float64)


It will help reader to always know what data type the array uses.

Please use the string version dtype='float64' (double quotes is OK if you like: dtype="float64") for fundamental types.

In terms of functionality it does not matter to use string or object for dtype. It's just for style consistency.

Implement LU factorization (PA=LU), linear system solve (Ax=b), and matrix inverse (A^-1) for float, double, complex64, and complex128 types. Expose lu_factorization, lu_solve, lu_inv as free functions. The .solve() and .inv() convenience methods on SimpleArray (for floating-point and complex types only) are attached from the linalg module via pybind11, so the buffer module has no dependency on linalg/. Tests cover factorization reconstruction (PA = LU) with explicit fixtures, known-solution solve/inverse checks, pivoting behavior, complex types, and error paths (non-square, 1D input, dimension mismatch, singular matrix). Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

tigercosmos · 2026-04-26T05:18:23Z

@yungyuc I have fixed the formating issue.

yungyuc

LGTM

tigercosmos force-pushed the linear_solver_0412 branch 3 times, most recently from ea0d54f to bc6585d Compare April 18, 2026 18:42

tigercosmos commented Apr 18, 2026

View reviewed changes

tigercosmos force-pushed the linear_solver_0412 branch from bc6585d to 8909a40 Compare April 18, 2026 19:01

tigercosmos force-pushed the linear_solver_0412 branch from 9f38505 to 862f097 Compare April 18, 2026 20:02

yungyuc requested changes Apr 19, 2026

View reviewed changes

yungyuc requested review from ThreeMonth03 and j8xixo12 April 19, 2026 01:10

yungyuc assigned tigercosmos Apr 19, 2026

yungyuc added the array Multi-dimensional array implementation label Apr 19, 2026

yungyuc added this to tensor operations Apr 19, 2026

github-project-automation Bot moved this to Todo in tensor operations Apr 19, 2026

yungyuc moved this from Todo to In Progress in tensor operations Apr 19, 2026

tigercosmos marked this pull request as draft April 19, 2026 10:24

tigercosmos force-pushed the linear_solver_0412 branch 4 times, most recently from 83762a9 to e758716 Compare April 19, 2026 12:19

tigercosmos marked this pull request as ready for review April 19, 2026 12:24

yungyuc requested changes Apr 20, 2026

View reviewed changes

tigercosmos force-pushed the linear_solver_0412 branch from e758716 to 7eba43d Compare April 22, 2026 18:18

tigercosmos force-pushed the linear_solver_0412 branch from 7eba43d to c8ced39 Compare April 22, 2026 18:29

tigercosmos commented Apr 22, 2026

View reviewed changes

yungyuc requested changes Apr 26, 2026

View reviewed changes

tigercosmos force-pushed the linear_solver_0412 branch from c8ced39 to 0eb0621 Compare April 26, 2026 05:18

yungyuc approved these changes Apr 26, 2026

View reviewed changes

yungyuc merged commit 3d73a5d into solvcon:master Apr 26, 2026
14 checks passed

github-project-automation Bot moved this from In Progress to Done in tensor operations Apr 26, 2026

tigercosmos deleted the linear_solver_0412 branch April 26, 2026 06:15

yungyuc mentioned this pull request Apr 27, 2026

Implement a linear solver for SimpleArray #719

Closed

yungyuc mentioned this pull request May 23, 2026

Make the LU factorization helper class modmesh::detail::Lu stateful and outside detail #809

Open

		np.testing.assert_allclose(ps_upd, ps_upd_np, atol=1e-12, rtol=0.0)


		class TestLuSolver(unittest.TestCase):

Conversation

tigercosmos commented Apr 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Test plan

Document

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tigercosmos commented Apr 18, 2026

Uh oh!

yungyuc left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yungyuc commented Apr 19, 2026

Document

Uh oh!

tigercosmos commented Apr 19, 2026

Uh oh!

yungyuc commented Apr 19, 2026

Uh oh!

tigercosmos commented Apr 19, 2026

Uh oh!

yungyuc commented Apr 19, 2026

Uh oh!

tigercosmos commented Apr 19, 2026

Uh oh!

tigercosmos commented Apr 19, 2026

Uh oh!

yungyuc left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tigercosmos left a comment

Choose a reason for hiding this comment

tigercosmos commented Apr 18, 2026 •

edited

Loading

yungyuc left a comment •

edited

Loading

yungyuc left a comment •

edited

Loading

yungyuc left a comment •

edited

Loading