shithub: amd64-simd

Download patch

ref: 39bdd5dc96197496ebb042fa1f1ad1a52738b33f
parent: 334e0b89bcc263c5155dfa9829aaaaed4c35d3f2
author: rodri <rgl@antares-labs.eu>
date: Sun Apr 6 11:23:08 EDT 2025 

bench: fix mulm3?_unrl and optimize harder


--- a/bench/main.c
+++ b/bench/main.c
@@ -96,29 +96,20 @@
 void
 mulm_unrl(Matrix a, Matrix b)
 {
-	Matrix tmp;
+	double t0, t1, t2;
 
-	tmp[0][0] = 0;
-	tmp[0][0] += a[0][0]*b[0][0]; tmp[0][0] += a[0][1]*b[1][0]; tmp[0][0] += a[0][2]*b[2][0];
-	tmp[0][1] = 0;
-	tmp[0][1] += a[0][0]*b[0][1]; tmp[0][1] += a[0][1]*b[1][1]; tmp[0][1] += a[0][2]*b[2][1];
-	tmp[0][2] = 0;
-	tmp[0][2] += a[0][0]*b[0][2]; tmp[0][2] += a[0][1]*b[1][2]; tmp[0][2] += a[0][2]*b[2][2];
-
-	tmp[1][0] = 0;
-	tmp[1][0] += a[1][0]*b[0][0]; tmp[1][0] += a[1][1]*b[1][0]; tmp[1][0] += a[1][2]*b[2][0];
-	tmp[1][1] = 0;
-	tmp[1][1] += a[1][0]*b[0][1]; tmp[1][1] += a[1][1]*b[1][1]; tmp[1][1] += a[1][2]*b[2][1];
-	tmp[1][2] = 0;
-	tmp[1][2] += a[1][0]*b[0][2]; tmp[1][2] += a[1][1]*b[1][2]; tmp[1][2] += a[1][2]*b[2][2];
-
-	tmp[2][0] = 0;
-	tmp[2][0] += a[2][0]*b[0][0]; tmp[2][0] += a[2][1]*b[1][0]; tmp[2][0] += a[2][2]*b[2][0];
-	tmp[2][1] = 0;
-	tmp[2][1] += a[2][0]*b[0][1]; tmp[2][1] += a[2][1]*b[1][1]; tmp[2][1] += a[2][2]*b[2][1];
-	tmp[2][2] = 0;
-	tmp[2][2] += a[2][0]*b[0][2]; tmp[2][2] += a[2][1]*b[1][2]; tmp[2][2] += a[2][2]*b[2][2];
-	memmove(a, tmp, 3*3*sizeof(double));
+	t0 = a[0][0]; t1 = a[0][1]; t2 = a[0][2];
+	a[0][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0];
+	a[0][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1];
+	a[0][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2];
+	t0 = a[1][0]; t1 = a[1][1]; t2 = a[1][2];
+	a[1][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0];
+	a[1][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1];
+	a[1][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2];
+	t0 = a[2][0]; t1 = a[2][1]; t2 = a[2][2];
+	a[2][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0];
+	a[2][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1];
+	a[2][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2];
 }
 
 void
@@ -162,60 +153,28 @@
 void
 mulm3_unrl(Matrix3 a, Matrix3 b)
 {
-	Matrix3 tmp;
+	double t0, t1, t2, t3;
 
-	tmp[0][0] = 0;
-	tmp[0][0] += a[0][0]*b[0][0]; tmp[0][0] += a[0][1]*b[1][0];
-	tmp[0][0] += a[0][2]*b[2][0]; tmp[0][0] += a[0][3]*b[3][0];
-	tmp[0][1] = 0;
-	tmp[0][1] += a[0][0]*b[0][1]; tmp[0][1] += a[0][1]*b[1][1];
-	tmp[0][1] += a[0][2]*b[2][1]; tmp[0][1] += a[0][3]*b[3][1];
-	tmp[0][2] = 0;
-	tmp[0][2] += a[0][0]*b[0][2]; tmp[0][2] += a[0][1]*b[1][2];
-	tmp[0][2] += a[0][2]*b[2][2]; tmp[0][2] += a[0][3]*b[3][2];
-	tmp[0][3] = 0;
-	tmp[0][3] += a[0][0]*b[0][3]; tmp[0][3] += a[0][1]*b[1][3];
-	tmp[0][3] += a[0][2]*b[2][3]; tmp[0][3] += a[0][3]*b[3][3];
-
-	tmp[1][0] = 0;
-	tmp[1][0] += a[1][0]*b[0][0]; tmp[1][0] += a[1][1]*b[1][0];
-	tmp[1][0] += a[1][2]*b[2][0]; tmp[1][0] += a[1][3]*b[3][0];
-	tmp[1][1] = 0;
-	tmp[1][1] += a[1][0]*b[0][1]; tmp[1][1] += a[1][1]*b[1][1];
-	tmp[1][1] += a[1][2]*b[2][1]; tmp[1][1] += a[1][3]*b[3][1];
-	tmp[1][2] = 0;
-	tmp[1][2] += a[1][0]*b[0][2]; tmp[1][2] += a[1][1]*b[1][2];
-	tmp[1][2] += a[1][2]*b[2][2]; tmp[1][2] += a[1][3]*b[3][2];
-	tmp[1][3] = 0;
-	tmp[1][3] += a[1][0]*b[0][3]; tmp[1][3] += a[1][1]*b[1][3];
-	tmp[1][3] += a[1][2]*b[2][3]; tmp[1][3] += a[1][3]*b[3][3];
-
-	tmp[2][0] = 0;
-	tmp[2][0] += a[2][0]*b[0][0]; tmp[2][0] += a[2][1]*b[1][0];
-	tmp[2][0] += a[2][2]*b[2][0]; tmp[2][0] += a[2][3]*b[3][0];
-	tmp[2][1] = 0;
-	tmp[2][1] += a[2][0]*b[0][1]; tmp[2][1] += a[2][1]*b[1][1];
-	tmp[2][1] += a[2][2]*b[2][1]; tmp[2][1] += a[2][3]*b[3][1];
-	tmp[2][2] = 0;
-	tmp[2][2] += a[2][0]*b[0][2]; tmp[2][2] += a[2][1]*b[1][2];
-	tmp[2][2] += a[2][2]*b[2][2]; tmp[2][2] += a[2][3]*b[3][2];
-	tmp[2][3] = 0;
-	tmp[2][3] += a[2][0]*b[0][3]; tmp[2][3] += a[2][1]*b[1][3];
-	tmp[2][3] += a[2][2]*b[2][3]; tmp[2][3] += a[2][3]*b[3][3];
-
-	tmp[3][0] = 0;
-	tmp[3][0] += a[3][0]*b[0][0]; tmp[3][0] += a[3][1]*b[1][0];
-	tmp[3][0] += a[3][2]*b[2][0]; tmp[3][0] += a[3][3]*b[3][0];
-	tmp[3][1] = 0;
-	tmp[3][1] += a[3][0]*b[0][1]; tmp[3][1] += a[3][1]*b[1][1];
-	tmp[3][1] += a[3][2]*b[2][1]; tmp[3][1] += a[3][3]*b[3][1];
-	tmp[3][2] = 0;
-	tmp[3][2] += a[3][0]*b[0][2]; tmp[3][2] += a[3][1]*b[1][2];
-	tmp[3][2] += a[3][2]*b[2][2]; tmp[3][2] += a[3][3]*b[3][2];
-	tmp[3][3] = 0;
-	tmp[3][3] += a[3][0]*b[0][3]; tmp[3][3] += a[3][1]*b[1][3];
-	tmp[3][3] += a[3][2]*b[2][3]; tmp[3][3] += a[3][3]*b[3][3];
-	memmove(a, tmp, 4*4*sizeof(double));
+	t0 = a[0][0]; t1 = a[0][1]; t2 = a[0][2]; t3 = a[0][3];
+	a[0][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0] + t3*b[3][0];
+	a[0][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1] + t3*b[3][1];
+	a[0][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2] + t3*b[3][2];
+	a[0][3] = t0*b[0][3] + t1*b[1][3] + t2*b[2][3] + t3*b[3][3];
+	t0 = a[1][0]; t1 = a[1][1]; t2 = a[1][2]; t3 = a[1][3];
+	a[1][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0] + t3*b[3][0];
+	a[1][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1] + t3*b[3][1];
+	a[1][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2] + t3*b[3][2];
+	a[1][3] = t0*b[0][3] + t1*b[1][3] + t2*b[2][3] + t3*b[3][3];
+	t0 = a[2][0]; t1 = a[2][1]; t2 = a[2][2]; t3 = a[2][3];
+	a[2][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0] + t3*b[3][0];
+	a[2][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1] + t3*b[3][1];
+	a[2][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2] + t3*b[3][2];
+	a[2][3] = t0*b[0][3] + t1*b[1][3] + t2*b[2][3] + t3*b[3][3];
+	t0 = a[3][0]; t1 = a[3][1]; t2 = a[3][2]; t3 = a[3][3];
+	a[3][0] = t0*b[0][0] + t1*b[1][0] + t2*b[2][0] + t3*b[3][0];
+	a[3][1] = t0*b[0][1] + t1*b[1][1] + t2*b[2][1] + t3*b[3][1];
+	a[3][2] = t0*b[0][2] + t1*b[1][2] + t2*b[2][2] + t3*b[3][2];
+	a[3][3] = t0*b[0][3] + t1*b[1][3] + t2*b[2][3] + t3*b[3][3];
 }
 
 static void
@@ -556,6 +515,30 @@
 	benchfreegr(&g);
 }
 
+static int
+eqmat2(Matrix a, Matrix b)
+{
+	int i, j;
+
+	for(i = 0; i < 3; i++)
+	for(j = 0; j < 3; j++)
+		if(a[i][j] != b[i][j])
+			return 0;
+	return 1;
+}
+
+static int
+eqmat3(Matrix3 a, Matrix3 b)
+{
+	int i, j;
+
+	for(i = 0; i < 4; i++)
+	for(j = 0; j < 4; j++)
+		if(a[i][j] != b[i][j])
+			return 0;
+	return 1;
+}
+
 static void
 bmulm(int fd)
 {
@@ -581,6 +564,14 @@
 		benchin(b0);
 		for(i = 0; i < 1e6; i++){
 			mulm(a, b);
+//			if(i == 0){
+//				Matrix t;
+//				memmove(t, a, 3*3*sizeof(double));
+//				memmove(a, a0, 3*3*sizeof(double));
+//				mulm_unrl(a, b);
+//				print("match %d\n", eqmat2(t, a));
+//				return;
+//			}
 			memmove(a, a0, 3*3*sizeof(double));
 		}
 		benchout(b0);
@@ -636,6 +627,14 @@
 		benchin(b0);
 		for(i = 0; i < 1e6; i++){
 			mulm3(a, b);
+//			if(i == 0){
+//				Matrix3 t;
+//				memmove(t, a, 4*4*sizeof(double));
+//				memmove(a, a0, 4*4*sizeof(double));
+//				mulm3_unrl(a, b);
+//				print("match %d\n", eqmat3(t, a));
+//				return;
+//			}
 			memmove(a, a0, 4*4*sizeof(double));
 		}
 		benchout(b0);
-- 
⑨