From e5ea80d0931fa3213ca9af3c9a203c14a3ce0689 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 5 Oct 2026 13:46:27 +0000 Subject: [PATCH 1/2] perf(math): build matrix results without zero-filling them first Every Matrix operator that returns a new matrix value-initialised its result (data = {}) and then overwrote every element. Without the O3 pragmas #348 removed, a consumer that builds at -O2 keeps those loops rolled, so the zeroing survives. For a 3x3 result it is a memset call, and the Arm GNU toolchain's newlib-nano memset stores one byte at a time. e-foc runs an RLS update in the control interrupt during mechanical identification. Its covariance update, (P - g g^T / d) / lambda, builds four 3x3 temporaries. In e-foc's SIL scenario that times the control interrupt through a full calibration, the slowest execution rose from 237 to 377 cycles of the emulated 25 MHz clock (arm-none-eabi-gcc 15.2.1). The four memset calls were 596 of its 1892 instructions. - Matrix gains a private constructor that leaves the elements uninitialised. operator+, operator-, the matrix and scalar products, Transpose, GetBlock and GetColumn build their results with it, since each writes every element. - The default and initializer-list constructors still zero-fill, so a default-constructed or partly listed matrix keeps its zeros. - A test evaluates those operators in constant expressions, so a compiler that checks (Clang does, GCC does not) rejects an operator that leaves an element unwritten. The RLS update keeps its matrix expression and now makes no memset calls. The same scenario peaks at 250 cycles. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HKhRVwu8yvxdKJgNo17oNc --- numerical/math/Matrix.hpp | 30 +++++++++++++++++++++--------- numerical/math/test/TestMatrix.cpp | 18 ++++++++++++++++++ 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/numerical/math/Matrix.hpp b/numerical/math/Matrix.hpp index 18bca0e..80fc7f9 100644 --- a/numerical/math/Matrix.hpp +++ b/numerical/math/Matrix.hpp @@ -63,7 +63,7 @@ namespace math [[nodiscard]] OPTIMIZE_FOR_SPEED friend constexpr Matrix operator+(const Matrix& lhs, const Matrix& rhs) { - Matrix result; + Matrix result{ Uninitialized{} }; for (size_type i = 0; i < size; ++i) result.data[i] = lhs.data[i] + rhs.data[i]; return result; @@ -71,7 +71,7 @@ namespace math [[nodiscard]] OPTIMIZE_FOR_SPEED friend constexpr Matrix operator-(const Matrix& lhs, const Matrix& rhs) { - Matrix result; + Matrix result{ Uninitialized{} }; for (size_type i = 0; i < size; ++i) result.data[i] = lhs.data[i] - rhs.data[i]; return result; @@ -80,7 +80,7 @@ namespace math template [[nodiscard]] OPTIMIZE_FOR_SPEED friend constexpr Matrix operator*(const Matrix& lhs, const Matrix& rhs) { - Matrix result; + Matrix result{ typename Matrix::Uninitialized{} }; for (size_type i = 0; i < Rows; ++i) { for (size_type j = 0; j < RhsCols; ++j) @@ -96,7 +96,7 @@ namespace math [[nodiscard]] OPTIMIZE_FOR_SPEED friend constexpr Matrix operator*(const Matrix& lhs, const T& scalar) { - Matrix result; + Matrix result{ Uninitialized{} }; for (size_type i = 0; i < size; ++i) result.data[i] = lhs.data[i] * scalar; return result; @@ -121,7 +121,16 @@ namespace math [[nodiscard]] constexpr Matrix GetColumn(size_type col) const; private: - std::array data = {}; + template + friend class Matrix; + + struct Uninitialized + {}; + + constexpr explicit Matrix(Uninitialized) noexcept + {} + + std::array data; }; template @@ -141,10 +150,13 @@ namespace math } template - constexpr Matrix::Matrix() noexcept = default; + constexpr Matrix::Matrix() noexcept + : data{} + {} template OPTIMIZE_FOR_SPEED constexpr Matrix::Matrix(std::initializer_list> init) + : data{} { size_t row = 0; for (const auto& row_list : init) @@ -267,7 +279,7 @@ namespace math OPTIMIZE_FOR_SPEED constexpr Matrix Matrix::Transpose() const { - Matrix result; + Matrix result{ typename Matrix::Uninitialized{} }; for (size_type i = 0; i < Rows; ++i) for (size_type j = 0; j < Cols; ++j) result.at(j, i) = at(i, j); @@ -328,7 +340,7 @@ namespace math { static_assert(BlockRows <= Rows && BlockCols <= Cols, "Requested block exceeds source matrix dimensions"); - Matrix result; + Matrix result{ typename Matrix::Uninitialized{} }; for (size_type r = 0; r < BlockRows; ++r) for (size_type c = 0; c < BlockCols; ++c) result.at(r, c) = at(rowOffset + r, colOffset + c); @@ -339,7 +351,7 @@ namespace math [[nodiscard]] OPTIMIZE_FOR_SPEED constexpr Matrix Matrix::GetColumn(size_type col) const { - Vector result; + Vector result{ typename Vector::Uninitialized{} }; for (size_type r = 0; r < Rows; ++r) result.at(r, 0) = at(r, col); return result; diff --git a/numerical/math/test/TestMatrix.cpp b/numerical/math/test/TestMatrix.cpp index 163a019..0b2ec28 100644 --- a/numerical/math/test/TestMatrix.cpp +++ b/numerical/math/test/TestMatrix.cpp @@ -427,6 +427,24 @@ TEST_F(MatrixFloatShapeTest, three_by_three_compound_assignment_and_trace) EXPECT_NEAR(square3.Transpose().at(0, 2), 7.0f, math::Tolerance()); } +TEST_F(MatrixFloatShapeTest, operators_write_every_element_of_their_result) +{ + constexpr math::Matrix column{ { 1.0f }, { 2.0f }, { 3.0f } }; + constexpr auto outer = column * column.Transpose(); + constexpr auto sum = outer + outer; + constexpr auto difference = sum - outer; + constexpr auto scaled = difference * 0.5f; + constexpr auto block = scaled.GetBlock<2, 2>(1, 1); + constexpr auto lastColumn = scaled.GetColumn(2); + + EXPECT_NEAR(outer.at(2, 1), 6.0f, math::Tolerance()); + EXPECT_NEAR(sum.at(0, 2), 6.0f, math::Tolerance()); + EXPECT_NEAR(difference.at(1, 1), 4.0f, math::Tolerance()); + EXPECT_NEAR(scaled.at(2, 2), 4.5f, math::Tolerance()); + EXPECT_NEAR(block.at(1, 0), 3.0f, math::Tolerance()); + EXPECT_NEAR(lastColumn.at(0, 0), 1.5f, math::Tolerance()); +} + TEST_F(MatrixFloatShapeTest, non_square_shapes_transpose_and_index) { math::Matrix row{ { 1.0f, 2.0f, 3.0f, 4.0f } }; From b54894adc0f1390c2f5ed340a4444fb3e91f9eb8 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 5 Oct 2026 14:34:17 +0000 Subject: [PATCH 2/2] fix(math): build the matrix product's result in a member function The matrix product is a friend function of its left operand's Matrix specialisation. It built its result, a Matrix of another shape, with that shape's private constructor. GCC and Clang allowed it, but MSVC rejects it (C2248): friendship goes to the Matrix specialisations, not to their friend functions. The operator now calls a private member function, Multiply, which builds the result in place, as Transpose, GetBlock and GetColumn already do. Members of every Matrix specialisation are friends of the others, so every compiler accepts it. A static helper that returned the uninitialised result was not used: GCC folds such an argument-less constexpr call into a zero-filled constant, which brings back the memset calls. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HKhRVwu8yvxdKJgNo17oNc --- numerical/math/Matrix.hpp | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/numerical/math/Matrix.hpp b/numerical/math/Matrix.hpp index 80fc7f9..9c10356 100644 --- a/numerical/math/Matrix.hpp +++ b/numerical/math/Matrix.hpp @@ -80,18 +80,7 @@ namespace math template [[nodiscard]] OPTIMIZE_FOR_SPEED friend constexpr Matrix operator*(const Matrix& lhs, const Matrix& rhs) { - Matrix result{ typename Matrix::Uninitialized{} }; - for (size_type i = 0; i < Rows; ++i) - { - for (size_type j = 0; j < RhsCols; ++j) - { - T sum{}; - for (size_type k = 0; k < Cols; ++k) - sum += lhs.at(i, k) * rhs.at(k, j); - result.at(i, j) = sum; - } - } - return result; + return lhs.Multiply(rhs); } [[nodiscard]] OPTIMIZE_FOR_SPEED friend constexpr Matrix operator*(const Matrix& lhs, const T& scalar) @@ -130,6 +119,9 @@ namespace math constexpr explicit Matrix(Uninitialized) noexcept {} + template + [[nodiscard]] constexpr Matrix Multiply(const Matrix& rhs) const; + std::array data; }; @@ -275,6 +267,25 @@ namespace math return *this; } + template + template + OPTIMIZE_FOR_SPEED constexpr Matrix + Matrix::Multiply(const Matrix& rhs) const + { + Matrix result{ typename Matrix::Uninitialized{} }; + for (size_type i = 0; i < Rows; ++i) + { + for (size_type j = 0; j < RhsCols; ++j) + { + T sum{}; + for (size_type k = 0; k < Cols; ++k) + sum += at(i, k) * rhs.at(k, j); + result.at(i, j) = sum; + } + } + return result; + } + template OPTIMIZE_FOR_SPEED constexpr Matrix Matrix::Transpose() const