From a574bec484a81e5c5bccaf18030d17056907a950 Mon Sep 17 00:00:00 2001
From: Frederick Roy <froy@lnrobo.com>
Date: Mon, 31 Aug 2026 10:24:24 +0900
Subject: [PATCH 1/3] demonstrate that addTransposed/subTransposed selfaliasing
 is wrong

---
 Sofa/framework/Type/test/MatTypes_test.cpp | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/Sofa/framework/Type/test/MatTypes_test.cpp b/Sofa/framework/Type/test/MatTypes_test.cpp
index 4251b0c7231..decf3caf89d 100644
--- a/Sofa/framework/Type/test/MatTypes_test.cpp
+++ b/Sofa/framework/Type/test/MatTypes_test.cpp
@@ -699,6 +699,28 @@ TEST(MatTypesTest, addSubTransposed)
     EXPECT_EQ(M, Matrix3::Identity() + A.transposed());
 }
 
+TEST(MatTypesTest, addTransposedSelfAliasing)
+{
+    Matrix3 M(Matrix3::Line(1., 2., 3.), Matrix3::Line(4., 5., 6.), Matrix3::Line(7., 8., 9.));
+    const Matrix3 expected(Matrix3::Line(2., 6., 10.),
+                           Matrix3::Line(6., 10., 14.),
+                           Matrix3::Line(10., 14., 18.));
+
+    M.addTransposed(M);
+    EXPECT_EQ(M, expected);
+}
+
+TEST(MatTypesTest, subTransposedSelfAliasing)
+{
+    Matrix3 M(Matrix3::Line(1., 2., 3.), Matrix3::Line(4., 5., 6.), Matrix3::Line(7., 8., 9.));
+    const Matrix3 expected(Matrix3::Line(0., -2., -4.),
+                           Matrix3::Line(2., 0., -2.),
+                           Matrix3::Line(4., 2., 0.));
+
+    M.subTransposed(M);
+    EXPECT_EQ(M, expected);
+}
+
 TEST(MatTypesTest, symmetrize)
 {
     Matrix3 A(Matrix3::Line(1., 2., 3.), Matrix3::Line(4., 5., 6.), Matrix3::Line(7., 8., 9.));

From 07fe9d4be511544d25d386baef8a5abf9b4e9dc3 Mon Sep 17 00:00:00 2001
From: Frederick Roy <froy@lnrobo.com>
Date: Mon, 31 Aug 2026 10:35:12 +0900
Subject: [PATCH 2/3] fixes selfaliases version

---
 Sofa/framework/Type/src/sofa/type/Mat.h | 49 ++++++++++++++++++++++---
 1 file changed, 43 insertions(+), 6 deletions(-)

diff --git a/Sofa/framework/Type/src/sofa/type/Mat.h b/Sofa/framework/Type/src/sofa/type/Mat.h
index 7f37337a5cf..0c8052882ac 100644
--- a/Sofa/framework/Type/src/sofa/type/Mat.h
+++ b/Sofa/framework/Type/src/sofa/type/Mat.h
@@ -708,17 +708,54 @@ class Mat
     /// Addition of the transposed of m
     constexpr void addTransposed(const Mat<C,L,real>& m) noexcept
     {
-        for(Size i=0; i<L; i++)
-            for(Size j=0; j<C; j++)
-                (*this)(i,j) += m(j,i);
+        if (canSelfTranspose(*this, m))
+        {
+            // m aliases *this: process each off-diagonal pair once, from values read
+            // before either of the two positions is written.
+            for(Size i=0; i<L; i++)
+            {
+                for(Size j=i+1; j<C; j++)
+                {
+                    const real sum = (*this)(i,j) + (*this)(j,i);
+                    (*this)(i,j) = sum;
+                    (*this)(j,i) = sum;
+                }
+                (*this)(i,i) += (*this)(i,i);
+            }
+        }
+        else
+        {
+            for(Size i=0; i<L; i++)
+                for(Size j=0; j<C; j++)
+                    (*this)(i,j) += m(j,i);
+        }
     }
 
     /// Subtraction of the transposed of m
     constexpr void subTransposed(const Mat<C,L,real>& m) noexcept
     {
-        for(Size i=0; i<L; i++)
-            for(Size j=0; j<C; j++)
-                (*this)(i,j) -= m(j,i);
+        if (canSelfTranspose(*this, m))
+        {
+            // m aliases *this: process each off-diagonal pair once, from values read
+            // before either of the two positions is written.
+            for(Size i=0; i<L; i++)
+            {
+                for(Size j=i+1; j<C; j++)
+                {
+                    const real mij = (*this)(i,j);
+                    const real mji = (*this)(j,i);
+                    (*this)(i,j) = mij - mji;
+                    (*this)(j,i) = mji - mij;
+                }
+                (*this)(i,i) = real{};
+            }
+        }
+        else
+        {
+            for(Size i=0; i<L; i++)
+                for(Size j=0; j<C; j++)
+                    (*this)(i,j) -= m(j,i);
+        }
     }
 
     /// Subtraction assignment operator.

From 4d2e85e36506f2fa2df7d3286853b4ffc44096b5 Mon Sep 17 00:00:00 2001
From: Frederick Roy <froy@lnrobo.com>
Date: Thu, 3 Sep 2026 07:48:33 +0900
Subject: [PATCH 3/3] apply suggestion

---
 Sofa/framework/Type/src/sofa/type/Mat.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Sofa/framework/Type/src/sofa/type/Mat.h b/Sofa/framework/Type/src/sofa/type/Mat.h
index 0c8052882ac..bfc52d5c6bf 100644
--- a/Sofa/framework/Type/src/sofa/type/Mat.h
+++ b/Sofa/framework/Type/src/sofa/type/Mat.h
@@ -714,13 +714,13 @@ class Mat
             // before either of the two positions is written.
             for(Size i=0; i<L; i++)
             {
+                (*this)(i,i) += (*this)(i,i);
                 for(Size j=i+1; j<C; j++)
                 {
                     const real sum = (*this)(i,j) + (*this)(j,i);
                     (*this)(i,j) = sum;
                     (*this)(j,i) = sum;
                 }
-                (*this)(i,i) += (*this)(i,i);
             }
         }
         else
@@ -740,6 +740,7 @@ class Mat
             // before either of the two positions is written.
             for(Size i=0; i<L; i++)
             {
+                (*this)(i,i) = real{};
                 for(Size j=i+1; j<C; j++)
                 {
                     const real mij = (*this)(i,j);
@@ -747,7 +748,6 @@ class Mat
                     (*this)(i,j) = mij - mji;
                     (*this)(j,i) = mji - mij;
                 }
-                (*this)(i,i) = real{};
             }
         }
         else