Integrate a matrix-algebra regression test.

1 files changed, 52 insertions, 18 deletions

diff --git a/test_matrix.c b/test_matrix.c
index bbbaffeb..0f8f17d1 100644
--- a/test_matrix.c
+++ b/test_matrix.c
@@ -8,6 +8,7 @@
 #include <stdbool.h>
 #include <stdlib.h>
 #include <stdio.h>
+#include <math.h>
 
 #include "matrix.h"
 
@@ -18,29 +19,62 @@
 #  define UNUSED
 #endif
 
-static double a[4][4] = {
-	{1, 2, 3, 4}, 
-	{5, 6, 7, 8},
-	{9, 10, 11, 12},
-	{13, 14, 15, 16}};
-static double ainv[4][4] = {
-    {-0.0028097262031538, 0.0073561782338441, 0.0040563991615249, -0.0086028511922152},
-    {-0.0018928965221056, 0.00016424848911584, 0.0070451831437969, -0.0053165351108071},
-    {0.031623219351165,	0.00010097243183351, -0.020259108723076, -0.011465083059922},
-    {0.031734962175727, -0.10943930681752, -0.52209748013199,1.5998018247738},
+static struct {
+    double mat[4][4];
+    double inv[4][4];
+} inverses[] = {
+    /* identity matrix is self-inverse */
+    {
+	.mat = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}},
+	.inv = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}},
+    },
+    /* inverse of a diagonal matrix has reciprocal values */
+    {
+	.mat = {{10,0,0,0}, {0,10,0,0}, {0,0,10,0}, {0,0,0,10}},
+	.inv = {{0.1,0,0,0}, {0,0.1,0,0}, {0,0,0.1,0}, {0,0,0,0.1}},
+    },
 };
 
-int main(int argc UNUSED, char *argv[] UNUSED)
+static void dump(const char *label, double m[4][4])
 {
-    double inv[4][4];
+    printf("%s:\n", label);
+    printf("%f, %f, %f, %f\n", m[0][0], m[0][1], m[0][2], m[0][3]);
+    printf("%f, %f, %f, %f\n", m[1][0], m[1][1], m[1][2], m[1][3]);
+    printf("%f, %f, %f, %f\n", m[2][0], m[2][1], m[2][2], m[2][3]);
+    printf("%f, %f, %f, %f\n", m[3][0], m[3][1], m[3][2], m[3][3]);
+}
+
+static bool check_diag(int n, double a[4][4], double b[4][4])
+{
+#define approx(x, y)	(fabs(x - y) < 0.0001)
 
-    matrix_invert(a, inv);
-    printf("%f %f %f %f\n",
-	   inv[0][0] - ainv[0][0],
-	   inv[1][1] - ainv[1][1],
-	   inv[2][2] - ainv[2][2],
-	   inv[3][3] - ainv[3][3]
+    if (approx(b[0][0], a[0][0]) && approx(b[1][1], a[1][1]) &&
+	approx(b[2][2], a[2][2]) && approx(b[3][3], a[3][3]))
+	return true;
+    
+    dump("a", a);
+    dump("b", b);
+    printf("Test %d residuals: %f %f %f %f\n",
+	   n,
+	   b[0][0] - a[0][0],
+	   b[1][1] - a[1][1],
+	   b[2][2] - a[2][2],
+	   b[3][3] - a[3][3]
 	);
+   return true;
+}
+
+int main(int argc UNUSED, char *argv[] UNUSED)
+{
+    size_t i;
+
+    for (i = 0; i < sizeof(inverses) / sizeof(inverses[0]); i++) {
+	double inverse[4][4];
+	matrix_invert(inverses[i].mat, inverse);
+	if (!check_diag(i, inverse, inverses[i].inv))
+	    break;
+    }
 
+    printf("Matrix-algebra regression test succeeded\n");
     exit(0);
 }