1 /*
2 * Unit test for matrix-algebra code
3 *
4 * Check examples computed at
5 * http://www.elektro-energetika.cz/calculations/matreg.php
6 *
7 * This file is Copyright (c) 2010 by the GPSD project
8 * SPDX-License-Identifier: BSD-2-clause
9 */
10 #include <stdlib.h>
11 #include <stdbool.h>
12 #include <stdlib.h>
13 #include <stdio.h>
14 #include <math.h>
15
16 #include "../compiler.h"
17 #include "../matrix.h"
18
19 static struct {
20 double mat[4][4];
21 double inv[4][4];
22 } inverses[] = {
23 /* identity matrix is self-inverse */
24 {
25 .mat = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}},
26 .inv = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}},
27 },
28 /* inverse of a diagonal matrix has reciprocal values */
29 {
30 .mat = {{10,0,0,0}, {0,10,0,0}, {0,0,10,0}, {0,0,0,10}},
31 .inv = {{0.1,0,0,0}, {0,0.1,0,0}, {0,0,0.1,0}, {0,0,0,0.1}},
32 },
33 /* random values with asymmetrical off-diagonal elements */
34 {
35 .mat = {{1,0,0,0}, {0,1,0,-2}, {0, 2,1,-4}, {0,0,0,1}},
36 .inv = {{1,0,0,0}, {0,1,0, 2}, {0,-2,1, 0}, {0,0,0,1}},
37 },
38 {
39 .mat = {{6,-4,1,-3},{-4,7,3,2},{1,3,6,-4},{-3,2,-4,6}},
40 .inv = {{14,34,-40,-31},{34,84,-99,-77},{-40,-99,117,91},{-31,-77,91,71}},
41 },
42 };
43
dump(const char * label,double m[4][4])44 static void dump(const char *label, double m[4][4])
45 {
46 printf("%s:\n", label);
47 printf("%f, %f, %f, %f\n", m[0][0], m[0][1], m[0][2], m[0][3]);
48 printf("%f, %f, %f, %f\n", m[1][0], m[1][1], m[1][2], m[1][3]);
49 printf("%f, %f, %f, %f\n", m[2][0], m[2][1], m[2][2], m[2][3]);
50 printf("%f, %f, %f, %f\n", m[3][0], m[3][1], m[3][2], m[3][3]);
51 }
52
check_diag(int n,double a[4][4],double b[4][4])53 static bool check_diag(int n, double a[4][4], double b[4][4])
54 {
55 #define approx(x, y) (fabs(x - y) < 0.0001)
56
57 if (approx(b[0][0], a[0][0]) && approx(b[1][1], a[1][1]) &&
58 approx(b[2][2], a[2][2]) && approx(b[3][3], a[3][3]))
59 return true;
60
61 dump("a", a);
62 dump("b", b);
63 printf("Test %d residuals: %f %f %f %f\n",
64 n,
65 b[0][0] - a[0][0],
66 b[1][1] - a[1][1],
67 b[2][2] - a[2][2],
68 b[3][3] - a[3][3]
69 );
70 return true;
71 }
72
main(int argc UNUSED,char * argv[]UNUSED)73 int main(int argc UNUSED, char *argv[] UNUSED)
74 {
75 unsigned int i;
76
77 for (i = 0; i < sizeof(inverses) / sizeof(inverses[0]); i++) {
78 double inverse[4][4];
79 if (!matrix_invert(inverses[i].mat, inverse)) {
80 printf("Vanishing determinant in test %u\n", i);
81 }
82 if (!check_diag(i, inverse, inverses[i].inv))
83 break;
84 }
85
86 printf("Matrix-algebra regression test succeeded\n");
87 exit(0);
88 }
89