From 87e4d97778560414420431ab8aa90540b70c9639 Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 10:39:47 +0100 Subject: [PATCH 1/7] Actually test eigenvalues in Test 1 --- tests/mat_4x4_eigenvalues.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/mat_4x4_eigenvalues.cpp b/tests/mat_4x4_eigenvalues.cpp index 698ca0a..13cd07f 100644 --- a/tests/mat_4x4_eigenvalues.cpp +++ b/tests/mat_4x4_eigenvalues.cpp @@ -27,7 +27,7 @@ int main() sm::vec, 4> eigenvalues = A.eigenvalues(); // Note: Due to numerical precision in polynomial root finding, - // we verify the characteristic polynomial is correct + // we verify the characteristic polynomial is correct. std::cout << " Computed Eigenvalues: "; for (size_t i = 0; i < 4; ++i) { std::cout << eigenvalues[i].real(); @@ -38,6 +38,15 @@ int main() } std::cout << "\n"; + if ((eigenvalues[0].real() - A[0]) > 1e-10 + || (eigenvalues[1].real() - A[5]) > 1e-10 + || (eigenvalues[2].real() - A[10]) > 1e-10 + || (eigenvalues[3].real() - A[15]) > 1e-10) { + std::cout << "Eigenvalues do not seem to be correct\n"; + --rtn; + } + + // For diagonal matrix, we verify trace and determinant are correct double trace = A.trace(); double det = A.determinant(); From 3c094f55d65edd3c8dd45810d695223d4a706304 Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 10:40:21 +0100 Subject: [PATCH 2/7] Removes extraneous code --- sm/mat.cppm | 53 ++++++++++++++++++----------------------------------- 1 file changed, 18 insertions(+), 35 deletions(-) diff --git a/sm/mat.cppm b/sm/mat.cppm index e8cedcf..5179731 100644 --- a/sm/mat.cppm +++ b/sm/mat.cppm @@ -869,8 +869,7 @@ export namespace sm for (std::uint32_t k = 1u; k <= Nr; ++k) { M = (*this) * M; // M_k = A * M_{k-1}, where M_0 = I - F trace = M.trace(); - coeffs[Nr - k] = -trace / F(k); + coeffs[Nr - k] = -(M.trace()) / F(k); if (k < Nr) { // M = M_k + c_{n-k} * I @@ -993,43 +992,27 @@ export namespace sm constexpr F_el eps = F_el{1e-14}; // needs to be F_el if F is complex sm::vec x = {}; - constexpr bool human_readable_version = false; - - if constexpr (human_readable_version) { - // Developed first with human-friendly matrix array indexing like A(r, c) - for (std::uint32_t i = (Nr - 1u); i != std::numeric_limits::max(); i--) { - if (std::abs ((*this)(i, i)) > eps) { // abs value of the triangular element - const std::uint32_t n = Nr - i - 1; // number of elements in the sum - F sum = F{0}; - for (std::uint32_t j = 0; j < n; ++j) { - sum += (*this)(i, i + n - j) * x[i + n - j]; - } - x[i] = ((*this)(i, Nc - 1) - sum ) / (*this)(i, i); - } // else x[i] remains 0. This could mean no solution if (*this)(i, Nc - 1) != 0 - } - } else { - // But compile this; Converting operator indexing to array indices as: (r, c) -> idx - // = (c * Nr) + r and removing need for the sum variable. - for (std::uint32_t i = (Nr - 1u); i != std::numeric_limits::max(); i--) { - const std::uint32_t de = (i * Nr) + i; // diagonal element - const std::uint32_t le = ((Nc - 1) * Nr) + i; // last element - if (std::abs (this->arr[de]) > eps) { - const std::uint32_t n = Nr - i - 1; - x[i] = this->arr[le]; - for (std::uint32_t j = 0; j < n; ++j) { - const std::uint32_t se = (Nr - 1 - j) * Nr + i; - x[i] -= this->arr[se] * x[i + n - j]; - } - x[i] /= this->arr[de]; + + for (std::uint32_t i = (Nr - 1u); i != std::numeric_limits::max(); i--) { + const std::uint32_t de = (i * Nr) + i; // diagonal element A(r, r) + const std::uint32_t le = ((Nc - 1) * Nr) + i; // last element A(r, Nc - 1) + if (std::abs (this->arr[de]) > eps) { // std::abs not sm::cem::abs as F may be complex + const std::uint32_t n = Nr - i - 1; + x[i] = this->arr[le]; + for (std::uint32_t j = 0; j < n; ++j) { + const std::uint32_t se = (Nr - 1 - j) * Nr + i; // sum element A(i, i + n - j) + x[i] -= this->arr[se] * x[i + n - j]; + } + x[i] /= this->arr[de]; + } else { + if constexpr (sm::is_complex::value) { + if (this->arr[le] != F{0}) { x[i] = F{std::numeric_limits::quiet_NaN()}; } } else { - if constexpr (sm::is_complex::value) { - if (this->arr[le] != F{0}) { x[i] = F{std::numeric_limits::quiet_NaN()}; } - } else { - if (this->arr[le] != F{0}) { x[i] = std::numeric_limits::quiet_NaN(); } - } + if (this->arr[le] != F{0}) { x[i] = std::numeric_limits::quiet_NaN(); } // else x[i] remains 0 } } } + return x; } From 139b5d08828abd32a1ed68b92e5e907d3c5b5d01 Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 10:43:35 +0100 Subject: [PATCH 3/7] back_substitution is const noexcept --- sm/mat.cppm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sm/mat.cppm b/sm/mat.cppm index 5179731..7a366fd 100644 --- a/sm/mat.cppm +++ b/sm/mat.cppm @@ -986,7 +986,7 @@ export namespace sm // Solve matrix by back-substitution. row_echelon_form is fine (don't need reduced row // echelon form). Nans in result mean there was no solution. template requires (Nc == Nr + 1) - sm::vec back_substitution() + sm::vec back_substitution() const noexcept { using F_el = typename value_type_of::type; constexpr F_el eps = F_el{1e-14}; // needs to be F_el if F is complex From e970d913290c62569c1b789b49d8b49e7710b663 Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 10:49:13 +0100 Subject: [PATCH 4/7] Some information for the future programmer --- sm/mat.cppm | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sm/mat.cppm b/sm/mat.cppm index 7a366fd..730bf45 100644 --- a/sm/mat.cppm +++ b/sm/mat.cppm @@ -1005,6 +1005,9 @@ export namespace sm } x[i] /= this->arr[de]; } else { + // arr[de] == 0. If the value in the right-most column is 0 AND the value in the + // diagonal is 0, then we have a free value and x[i] remains 0. If the value in + // the right most column is != 0, but arr[de] == 0, we have no solution. if constexpr (sm::is_complex::value) { if (this->arr[le] != F{0}) { x[i] = F{std::numeric_limits::quiet_NaN()}; } } else { From ae605cc97ebe76bba761be5d0ed82fb4e8f298e8 Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 11:01:06 +0100 Subject: [PATCH 5/7] Mostly just code comments, but also remove an extraneous line --- sm/mat.cppm | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/sm/mat.cppm b/sm/mat.cppm index 730bf45..8a2718b 100644 --- a/sm/mat.cppm +++ b/sm/mat.cppm @@ -983,8 +983,15 @@ export namespace sm return rref; } - // Solve matrix by back-substitution. row_echelon_form is fine (don't need reduced row - // echelon form). Nans in result mean there was no solution. + /*! + * Using back-substitution, solve the equation Ax = B, where A and B are represented in + * *this in the usual way, with A being a square matrix on the left and B the right-most + * column. *this should thus be the 'augmented matrix'. This matrix should have been reduced + * to row_echelon_form before calling back_substitution(). + * + * \return A vec of the solutions, x. NaNs in result mean there was no solution. Currently, + * I don't represent a free variable. Is that possible? I *could* use limits::max? + */ template requires (Nc == Nr + 1) sm::vec back_substitution() const noexcept { @@ -1000,8 +1007,7 @@ export namespace sm const std::uint32_t n = Nr - i - 1; x[i] = this->arr[le]; for (std::uint32_t j = 0; j < n; ++j) { - const std::uint32_t se = (Nr - 1 - j) * Nr + i; // sum element A(i, i + n - j) - x[i] -= this->arr[se] * x[i + n - j]; + x[i] -= this->arr[(Nr - 1 - j) * Nr + i] * x[i + n - j]; // sum element is A(i, i + n - j) } x[i] /= this->arr[de]; } else { From e84eec0e57a8709517e47562093a27d371082caa Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 11:09:10 +0100 Subject: [PATCH 6/7] No need to introduce a change in that line for this PR --- tests/mat_4x4_eigenvalues.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/mat_4x4_eigenvalues.cpp b/tests/mat_4x4_eigenvalues.cpp index 13cd07f..92807bc 100644 --- a/tests/mat_4x4_eigenvalues.cpp +++ b/tests/mat_4x4_eigenvalues.cpp @@ -27,7 +27,7 @@ int main() sm::vec, 4> eigenvalues = A.eigenvalues(); // Note: Due to numerical precision in polynomial root finding, - // we verify the characteristic polynomial is correct. + // we verify the characteristic polynomial is correct std::cout << " Computed Eigenvalues: "; for (size_t i = 0; i < 4; ++i) { std::cout << eigenvalues[i].real(); From 7a9292ad59662e8761cc773df512d86b8817954f Mon Sep 17 00:00:00 2001 From: Seb James Date: Sat, 16 May 2026 11:09:22 +0100 Subject: [PATCH 7/7] Prefere static_cast in this code --- sm/mat.cppm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sm/mat.cppm b/sm/mat.cppm index 8a2718b..b0e20b8 100644 --- a/sm/mat.cppm +++ b/sm/mat.cppm @@ -869,7 +869,7 @@ export namespace sm for (std::uint32_t k = 1u; k <= Nr; ++k) { M = (*this) * M; // M_k = A * M_{k-1}, where M_0 = I - coeffs[Nr - k] = -(M.trace()) / F(k); + coeffs[Nr - k] = -M.trace() / static_cast(k); if (k < Nr) { // M = M_k + c_{n-k} * I