summaryrefslogtreecommitdiff
path: root/tests/matrix-test.c
diff options
context:
space:
mode:
Diffstat (limited to 'tests/matrix-test.c')
-rw-r--r--tests/matrix-test.c419
1 files changed, 419 insertions, 0 deletions
diff --git a/tests/matrix-test.c b/tests/matrix-test.c
new file mode 100644
index 00000000..5b0513f3
--- /dev/null
+++ b/tests/matrix-test.c
@@ -0,0 +1,419 @@
+/*
+ * Copyright © 2012 Collabora, Ltd.
+ *
+ * Permission to use, copy, modify, distribute, and sell this software and
+ * its documentation for any purpose is hereby granted without fee, provided
+ * that the above copyright notice appear in all copies and that both that
+ * copyright notice and this permission notice appear in supporting
+ * documentation, and that the name of the copyright holders not be used in
+ * advertising or publicity pertaining to distribution of the software
+ * without specific, written prior permission. The copyright holders make
+ * no representations about the suitability of this software for any
+ * purpose. It is provided "as is" without express or implied warranty.
+ *
+ * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS
+ * SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS, IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
+ * RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF
+ * CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
+ * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <unistd.h>
+#include <signal.h>
+#include <time.h>
+
+#include "../shared/matrix.h"
+
+struct inverse_matrix {
+ double LU[16]; /* column-major */
+ unsigned perm[4]; /* permutation */
+};
+
+static struct timespec begin_time;
+
+static void
+reset_timer(void)
+{
+ clock_gettime(CLOCK_MONOTONIC, &begin_time);
+}
+
+static double
+read_timer(void)
+{
+ struct timespec t;
+
+ clock_gettime(CLOCK_MONOTONIC, &t);
+ return (double)(t.tv_sec - begin_time.tv_sec) +
+ 1e-9 * (t.tv_nsec - begin_time.tv_nsec);
+}
+
+static double
+det3x3(const float *c0, const float *c1, const float *c2)
+{
+ return (double)
+ c0[0] * c1[1] * c2[2] +
+ c1[0] * c2[1] * c0[2] +
+ c2[0] * c0[1] * c1[2] -
+ c0[2] * c1[1] * c2[0] -
+ c1[2] * c2[1] * c0[0] -
+ c2[2] * c0[1] * c1[0];
+}
+
+static double
+determinant(const struct weston_matrix *m)
+{
+ double det = 0;
+#if 1
+ /* develop on last row */
+ det -= m->d[3 + 0 * 4] * det3x3(&m->d[4], &m->d[8], &m->d[12]);
+ det += m->d[3 + 1 * 4] * det3x3(&m->d[0], &m->d[8], &m->d[12]);
+ det -= m->d[3 + 2 * 4] * det3x3(&m->d[0], &m->d[4], &m->d[12]);
+ det += m->d[3 + 3 * 4] * det3x3(&m->d[0], &m->d[4], &m->d[8]);
+#else
+ /* develop on first row */
+ det += m->d[0 + 0 * 4] * det3x3(&m->d[5], &m->d[9], &m->d[13]);
+ det -= m->d[0 + 1 * 4] * det3x3(&m->d[1], &m->d[9], &m->d[13]);
+ det += m->d[0 + 2 * 4] * det3x3(&m->d[1], &m->d[5], &m->d[13]);
+ det -= m->d[0 + 3 * 4] * det3x3(&m->d[1], &m->d[5], &m->d[9]);
+#endif
+ return det;
+}
+
+static void
+print_permutation_matrix(const struct inverse_matrix *m)
+{
+ const unsigned *p = m->perm;
+ const char *row[4] = {
+ "1 0 0 0\n",
+ "0 1 0 0\n",
+ "0 0 1 0\n",
+ "0 0 0 1\n"
+ };
+
+ printf(" P =\n%s%s%s%s", row[p[0]], row[p[1]], row[p[2]], row[p[3]]);
+}
+
+static void
+print_LU_decomposition(const struct inverse_matrix *m)
+{
+ unsigned r, c;
+
+ printf(" L "
+ " U\n");
+ for (r = 0; r < 4; ++r) {
+ double v;
+
+ for (c = 0; c < 4; ++c) {
+ if (c < r)
+ v = m->LU[r + c * 4];
+ else if (c == r)
+ v = 1.0;
+ else
+ v = 0.0;
+ printf(" %12.6f", v);
+ }
+
+ printf(" | ");
+
+ for (c = 0; c < 4; ++c) {
+ if (c >= r)
+ v = m->LU[r + c * 4];
+ else
+ v = 0.0;
+ printf(" %12.6f", v);
+ }
+ printf("\n");
+ }
+}
+
+static void
+print_inverse_data_matrix(const struct inverse_matrix *m)
+{
+ unsigned r, c;
+
+ for (r = 0; r < 4; ++r) {
+ for (c = 0; c < 4; ++c)
+ printf(" %12.6f", m->LU[r + c * 4]);
+ printf("\n");
+ }
+
+ printf("permutation: ");
+ for (r = 0; r < 4; ++r)
+ printf(" %u", m->perm[r]);
+ printf("\n");
+}
+
+static void
+print_matrix(const struct weston_matrix *m)
+{
+ unsigned r, c;
+
+ for (r = 0; r < 4; ++r) {
+ for (c = 0; c < 4; ++c)
+ printf(" %14.6e", m->d[r + c * 4]);
+ printf("\n");
+ }
+}
+
+static double
+frand(void)
+{
+ double r = random();
+ return r / (double)(RAND_MAX / 2) - 1.0f;
+}
+
+static void
+randomize_matrix(struct weston_matrix *m)
+{
+ unsigned i;
+ for (i = 0; i < 16; ++i)
+#if 1
+ m->d[i] = frand() * exp(10.0 * frand());
+#else
+ m->d[i] = frand();
+#endif
+}
+
+/* Take a matrix, compute inverse, multiply together
+ * and subtract the identity matrix to get the error matrix.
+ * Return the largest absolute value from the error matrix.
+ */
+static double
+test_inverse(struct weston_matrix *m)
+{
+ unsigned i;
+ struct inverse_matrix q;
+ double errsup = 0.0;
+
+ if (matrix_invert(q.LU, q.perm, m) != 0)
+ return INFINITY;
+
+ for (i = 0; i < 4; ++i)
+ inverse_transform(q.LU, q.perm, &m->d[i * 4]);
+
+ m->d[0] -= 1.0f;
+ m->d[5] -= 1.0f;
+ m->d[10] -= 1.0f;
+ m->d[15] -= 1.0f;
+
+ for (i = 0; i < 16; ++i) {
+ double err = fabs(m->d[i]);
+ if (err > errsup)
+ errsup = err;
+ }
+
+ return errsup;
+}
+
+enum {
+ TEST_OK,
+ TEST_NOT_INVERTIBLE_OK,
+ TEST_FAIL,
+ TEST_COUNT
+};
+
+static int
+test(void)
+{
+ struct weston_matrix m;
+ double det, errsup;
+
+ randomize_matrix(&m);
+ det = determinant(&m);
+
+ errsup = test_inverse(&m);
+ if (errsup < 1e-6)
+ return TEST_OK;
+
+ if (fabs(det) < 1e-5 && isinf(errsup))
+ return TEST_NOT_INVERTIBLE_OK;
+
+ printf("test fail, det: %g, error sup: %g\n", det, errsup);
+
+ return TEST_FAIL;
+}
+
+static int running;
+static void
+stopme(int n)
+{
+ running = 0;
+}
+
+static void
+test_loop_precision(void)
+{
+ int counts[TEST_COUNT] = { 0 };
+
+ printf("\nRunning a test loop for 10 seconds...\n");
+ running = 1;
+ alarm(10);
+ while (running) {
+ counts[test()]++;
+ }
+
+ printf("tests: %d ok, %d not invertible but ok, %d failed.\n"
+ "Total: %d iterations.\n",
+ counts[TEST_OK], counts[TEST_NOT_INVERTIBLE_OK],
+ counts[TEST_FAIL],
+ counts[TEST_OK] + counts[TEST_NOT_INVERTIBLE_OK] +
+ counts[TEST_FAIL]);
+}
+
+static void __attribute__((noinline))
+test_loop_speed_matrixvector(void)
+{
+ struct weston_matrix m;
+ struct weston_vector v = { { 0.5, 0.5, 0.5, 1.0 } };
+ unsigned long count = 0;
+ double t;
+
+ printf("\nRunning 3 s test on weston_matrix_transform()...\n");
+
+ weston_matrix_init(&m);
+
+ running = 1;
+ alarm(3);
+ reset_timer();
+ while (running) {
+ weston_matrix_transform(&m, &v);
+ count++;
+ }
+ t = read_timer();
+
+ printf("%lu iterations in %f seconds, avg. %.1f us/iter.\n",
+ count, t, 1e9 * t / count);
+}
+
+static void __attribute__((noinline))
+test_loop_speed_inversetransform(void)
+{
+ struct weston_matrix m;
+ struct inverse_matrix inv;
+ struct weston_vector v = { { 0.5, 0.5, 0.5, 1.0 } };
+ unsigned long count = 0;
+ double t;
+
+ printf("\nRunning 3 s test on inverse_transform()...\n");
+
+ weston_matrix_init(&m);
+ matrix_invert(inv.LU, inv.perm, &m);
+
+ running = 1;
+ alarm(3);
+ reset_timer();
+ while (running) {
+ inverse_transform(inv.LU, inv.perm, v.f);
+ count++;
+ }
+ t = read_timer();
+
+ printf("%lu iterations in %f seconds, avg. %.1f us/iter.\n",
+ count, t, 1e9 * t / count);
+}
+
+static void __attribute__((noinline))
+test_loop_speed_invert(void)
+{
+ struct weston_matrix m;
+ struct inverse_matrix inv;
+ unsigned long count = 0;
+ double t;
+
+ printf("\nRunning 3 s test on matrix_invert()...\n");
+
+ weston_matrix_init(&m);
+
+ running = 1;
+ alarm(3);
+ reset_timer();
+ while (running) {
+ matrix_invert(inv.LU, inv.perm, &m);
+ count++;
+ }
+ t = read_timer();
+
+ printf("%lu iterations in %f seconds, avg. %.1f ns/iter.\n",
+ count, t, 1e9 * t / count);
+}
+
+static void __attribute__((noinline))
+test_loop_speed_invert_explicit(void)
+{
+ struct weston_matrix m;
+ unsigned long count = 0;
+ double t;
+
+ printf("\nRunning 3 s test on weston_matrix_invert()...\n");
+
+ weston_matrix_init(&m);
+
+ running = 1;
+ alarm(3);
+ reset_timer();
+ while (running) {
+ weston_matrix_invert(&m, &m);
+ count++;
+ }
+ t = read_timer();
+
+ printf("%lu iterations in %f seconds, avg. %.1f ns/iter.\n",
+ count, t, 1e9 * t / count);
+}
+
+int main(void)
+{
+ struct sigaction ding;
+ struct weston_matrix M;
+ struct inverse_matrix Q;
+ int ret;
+ double errsup;
+ double det;
+
+ ding.sa_handler = stopme;
+ sigemptyset(&ding.sa_mask);
+ ding.sa_flags = 0;
+ sigaction(SIGALRM, &ding, NULL);
+
+ srandom(13);
+
+ M.d[0] = 3.0; M.d[4] = 17.0; M.d[8] = 10.0; M.d[12] = 0.0;
+ M.d[1] = 2.0; M.d[5] = 4.0; M.d[9] = -2.0; M.d[13] = 0.0;
+ M.d[2] = 6.0; M.d[6] = 18.0; M.d[10] = -12; M.d[14] = 0.0;
+ M.d[3] = 0.0; M.d[7] = 0.0; M.d[11] = 0.0; M.d[15] = 1.0;
+
+ ret = matrix_invert(Q.LU, Q.perm, &M);
+ printf("ret = %d\n", ret);
+ printf("det = %g\n\n", determinant(&M));
+
+ if (ret != 0)
+ return 1;
+
+ print_inverse_data_matrix(&Q);
+ printf("P * A = L * U\n");
+ print_permutation_matrix(&Q);
+ print_LU_decomposition(&Q);
+
+
+ printf("a random matrix:\n");
+ randomize_matrix(&M);
+ det = determinant(&M);
+ print_matrix(&M);
+ errsup = test_inverse(&M);
+ printf("\nThe matrix multiplied by its inverse, error:\n");
+ print_matrix(&M);
+ printf("max abs error: %g, original determinant %g\n", errsup, det);
+
+ test_loop_precision();
+ test_loop_speed_matrixvector();
+ test_loop_speed_inversetransform();
+ test_loop_speed_invert();
+ test_loop_speed_invert_explicit();
+
+ return 0;
+}
View Lorry Controller Status