specialize mat4 multiply routines

mul_m4_m4m4(R, A, B) gives us R = AB in general. Existing code assumed the worst, that A and B both alias the output R. For safety it makes internal copies of A and B before calculating & writing R.

This is the least common case. Usually all 3 matrices differ. Often we see M = AM or M = MB, but never M = MM.

With this revision mul_m4_m4m4 is called in exactly the same way but copies inputs only when needed. If you know the inputs are independent of the output use the "uniq" variant to skip the saftety checks.

2 files changed, 48 insertions, 23 deletions

diff --git a/source/blender/blenlib/BLI_math_matrix.h b/source/blender/blenlib/BLI_math_matrix.h
index 6388105032a..a5868e00c8e 100644
--- a/source/blender/blenlib/BLI_math_matrix.h
+++ b/source/blender/blenlib/BLI_math_matrix.h
@@ -72,6 +72,11 @@ void mul_m4_m4m3(float R[4][4], const float A[4][4], const float B[3][3]);
 void mul_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4]);
 void mul_m3_m3m4(float R[3][3], const float A[4][4], const float B[3][3]);
 
+/* special matrix multiplies */
+void mul_m4_m4m4_uniq(float R[4][4], const float A[4][4], const float B[4][4]); /* R <-- AB, R is neither A nor B */
+void mul_m4_m4_pre(float R[4][4], const float A[4][4]); /* R <-- AR */
+void mul_m4_m4_post(float R[4][4], const float B[4][4]); /* R <-- RB */
+
 /* mul_m3_series */
 void _va_mul_m3_series_3(float R[3][3], float M1[3][3], float M2[3][3]) ATTR_NONNULL();
 void _va_mul_m3_series_4(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3]) ATTR_NONNULL();
diff --git a/source/blender/blenlib/intern/math_matrix.c b/source/blender/blenlib/intern/math_matrix.c
index 87b757ef4df..b22061b4ea8 100644
--- a/source/blender/blenlib/intern/math_matrix.c
+++ b/source/blender/blenlib/intern/math_matrix.c
@@ -177,36 +177,56 @@ void swap_m4m4(float m1[4][4], float m2[4][4])
 
 /******************************** Arithmetic *********************************/
 
-void mul_m4_m4m4(float m1[4][4], const float m3_[4][4], const float m2_[4][4])
+void mul_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4])
 {
-	float m2[4][4], m3[4][4];
+	if (A == R)
+		mul_m4_m4_post(R, B);
+	else if (B == R)
+		mul_m4_m4_pre(R, A);
+	else
+		mul_m4_m4m4_uniq(R, A, B);
+}
 
-	/* copy so it works when m1 is the same pointer as m2 or m3 */
-	/* TODO: avoid copying when matrices are different */
-	copy_m4_m4(m2, m2_);
-	copy_m4_m4(m3, m3_);
+void mul_m4_m4m4_uniq(float R[4][4], const float A[4][4], const float B[4][4])
+{
+	BLI_assert(R != A && R != B);
 
-	/* matrix product: m1[j][k] = m2[j][i].m3[i][k] */
-	m1[0][0] = m2[0][0] * m3[0][0] + m2[0][1] * m3[1][0] + m2[0][2] * m3[2][0] + m2[0][3] * m3[3][0];
-	m1[0][1] = m2[0][0] * m3[0][1] + m2[0][1] * m3[1][1] + m2[0][2] * m3[2][1] + m2[0][3] * m3[3][1];
-	m1[0][2] = m2[0][0] * m3[0][2] + m2[0][1] * m3[1][2] + m2[0][2] * m3[2][2] + m2[0][3] * m3[3][2];
-	m1[0][3] = m2[0][0] * m3[0][3] + m2[0][1] * m3[1][3] + m2[0][2] * m3[2][3] + m2[0][3] * m3[3][3];
+	/* matrix product: R[j][k] = A[j][i] . B[i][k] */
+	R[0][0] = B[0][0] * A[0][0] + B[0][1] * A[1][0] + B[0][2] * A[2][0] + B[0][3] * A[3][0];
+	R[0][1] = B[0][0] * A[0][1] + B[0][1] * A[1][1] + B[0][2] * A[2][1] + B[0][3] * A[3][1];
+	R[0][2] = B[0][0] * A[0][2] + B[0][1] * A[1][2] + B[0][2] * A[2][2] + B[0][3] * A[3][2];
+	R[0][3] = B[0][0] * A[0][3] + B[0][1] * A[1][3] + B[0][2] * A[2][3] + B[0][3] * A[3][3];
 
-	m1[1][0] = m2[1][0] * m3[0][0] + m2[1][1] * m3[1][0] + m2[1][2] * m3[2][0] + m2[1][3] * m3[3][0];
-	m1[1][1] = m2[1][0] * m3[0][1] + m2[1][1] * m3[1][1] + m2[1][2] * m3[2][1] + m2[1][3] * m3[3][1];
-	m1[1][2] = m2[1][0] * m3[0][2] + m2[1][1] * m3[1][2] + m2[1][2] * m3[2][2] + m2[1][3] * m3[3][2];
-	m1[1][3] = m2[1][0] * m3[0][3] + m2[1][1] * m3[1][3] + m2[1][2] * m3[2][3] + m2[1][3] * m3[3][3];
+	R[1][0] = B[1][0] * A[0][0] + B[1][1] * A[1][0] + B[1][2] * A[2][0] + B[1][3] * A[3][0];
+	R[1][1] = B[1][0] * A[0][1] + B[1][1] * A[1][1] + B[1][2] * A[2][1] + B[1][3] * A[3][1];
+	R[1][2] = B[1][0] * A[0][2] + B[1][1] * A[1][2] + B[1][2] * A[2][2] + B[1][3] * A[3][2];
+	R[1][3] = B[1][0] * A[0][3] + B[1][1] * A[1][3] + B[1][2] * A[2][3] + B[1][3] * A[3][3];
 
-	m1[2][0] = m2[2][0] * m3[0][0] + m2[2][1] * m3[1][0] + m2[2][2] * m3[2][0] + m2[2][3] * m3[3][0];
-	m1[2][1] = m2[2][0] * m3[0][1] + m2[2][1] * m3[1][1] + m2[2][2] * m3[2][1] + m2[2][3] * m3[3][1];
-	m1[2][2] = m2[2][0] * m3[0][2] + m2[2][1] * m3[1][2] + m2[2][2] * m3[2][2] + m2[2][3] * m3[3][2];
-	m1[2][3] = m2[2][0] * m3[0][3] + m2[2][1] * m3[1][3] + m2[2][2] * m3[2][3] + m2[2][3] * m3[3][3];
+	R[2][0] = B[2][0] * A[0][0] + B[2][1] * A[1][0] + B[2][2] * A[2][0] + B[2][3] * A[3][0];
+	R[2][1] = B[2][0] * A[0][1] + B[2][1] * A[1][1] + B[2][2] * A[2][1] + B[2][3] * A[3][1];
+	R[2][2] = B[2][0] * A[0][2] + B[2][1] * A[1][2] + B[2][2] * A[2][2] + B[2][3] * A[3][2];
+	R[2][3] = B[2][0] * A[0][3] + B[2][1] * A[1][3] + B[2][2] * A[2][3] + B[2][3] * A[3][3];
 
-	m1[3][0] = m2[3][0] * m3[0][0] + m2[3][1] * m3[1][0] + m2[3][2] * m3[2][0] + m2[3][3] * m3[3][0];
-	m1[3][1] = m2[3][0] * m3[0][1] + m2[3][1] * m3[1][1] + m2[3][2] * m3[2][1] + m2[3][3] * m3[3][1];
-	m1[3][2] = m2[3][0] * m3[0][2] + m2[3][1] * m3[1][2] + m2[3][2] * m3[2][2] + m2[3][3] * m3[3][2];
-	m1[3][3] = m2[3][0] * m3[0][3] + m2[3][1] * m3[1][3] + m2[3][2] * m3[2][3] + m2[3][3] * m3[3][3];
+	R[3][0] = B[3][0] * A[0][0] + B[3][1] * A[1][0] + B[3][2] * A[2][0] + B[3][3] * A[3][0];
+	R[3][1] = B[3][0] * A[0][1] + B[3][1] * A[1][1] + B[3][2] * A[2][1] + B[3][3] * A[3][1];
+	R[3][2] = B[3][0] * A[0][2] + B[3][1] * A[1][2] + B[3][2] * A[2][2] + B[3][3] * A[3][2];
+	R[3][3] = B[3][0] * A[0][3] + B[3][1] * A[1][3] + B[3][2] * A[2][3] + B[3][3] * A[3][3];
+}
+
+void mul_m4_m4_pre(float R[4][4], const float A[4][4])
+{
+	BLI_assert(A != R);
+	float B[4][4];
+	copy_m4_m4(B, R);
+	mul_m4_m4m4_uniq(R, A, B);
+}
 
+void mul_m4_m4_post(float R[4][4], const float B[4][4])
+{
+	BLI_assert(B != R);
+	float A[4][4];
+	copy_m4_m4(A, R);
+	mul_m4_m4m4_uniq(R, A, B);
 }
 
 void mul_m3_m3m3(float m1[3][3], const float m3_[3][3], const float m2_[3][3])