tests/matrix-test.c

   1 /*
   2  * Copyright © 2012 Collabora, Ltd.
   3  *
   4  * Permission to use, copy, modify, distribute, and sell this software and
   5  * its documentation for any purpose is hereby granted without fee, provided
   6  * that the above copyright notice appear in all copies and that both that
   7  * copyright notice and this permission notice appear in supporting
   8  * documentation, and that the name of the copyright holders not be used in
   9  * advertising or publicity pertaining to distribution of the software
  10  * without specific, written prior permission.  The copyright holders make
  11  * no representations about the suitability of this software for any
  12  * purpose.  It is provided "as is" without express or implied warranty.
  13  *
  14  * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS
  15  * SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
  16  * FITNESS, IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
  17  * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
  18  * RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF
  19  * CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
  20  * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  21  */
  22
  23 #include "config.h"
  24
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27 #include <math.h>
  28 #include <unistd.h>
  29 #include <signal.h>
  30 #include <time.h>
  31
  32 #include "../shared/matrix.h"
  33
  34 struct inverse_matrix {
  35         double LU[16];          /* column-major */
  36         unsigned perm[4];       /* permutation */
  37 };
  38
  39 static struct timespec begin_time;
  40
  41 static void
  42 reset_timer(void)
  43 {
  44         clock_gettime(CLOCK_MONOTONIC, &begin_time);
  45 }
  46
  47 static double
  48 read_timer(void)
  49 {
  50         struct timespec t;
  51
  52         clock_gettime(CLOCK_MONOTONIC, &t);
  53         return (double)(t.tv_sec - begin_time.tv_sec) +
  54                1e-9 * (t.tv_nsec - begin_time.tv_nsec);
  55 }
  56
  57 static double
  58 det3x3(const float *c0, const float *c1, const float *c2)
  59 {
  60         return (double)
  61                 c0[0] * c1[1] * c2[2] +
  62                 c1[0] * c2[1] * c0[2] +
  63                 c2[0] * c0[1] * c1[2] -
  64                 c0[2] * c1[1] * c2[0] -
  65                 c1[2] * c2[1] * c0[0] -
  66                 c2[2] * c0[1] * c1[0];
  67 }
  68
  69 static double
  70 determinant(const struct weston_matrix *m)
  71 {
  72         double det = 0;
  73 #if 1
  74         /* develop on last row */
  75         det -= m->d[3 + 0 * 4] * det3x3(&m->d[4], &m->d[8], &m->d[12]);
  76         det += m->d[3 + 1 * 4] * det3x3(&m->d[0], &m->d[8], &m->d[12]);
  77         det -= m->d[3 + 2 * 4] * det3x3(&m->d[0], &m->d[4], &m->d[12]);
  78         det += m->d[3 + 3 * 4] * det3x3(&m->d[0], &m->d[4], &m->d[8]);
  79 #else
  80         /* develop on first row */
  81         det += m->d[0 + 0 * 4] * det3x3(&m->d[5], &m->d[9], &m->d[13]);
  82         det -= m->d[0 + 1 * 4] * det3x3(&m->d[1], &m->d[9], &m->d[13]);
  83         det += m->d[0 + 2 * 4] * det3x3(&m->d[1], &m->d[5], &m->d[13]);
  84         det -= m->d[0 + 3 * 4] * det3x3(&m->d[1], &m->d[5], &m->d[9]);
  85 #endif
  86         return det;
  87 }
  88
  89 static void
  90 print_permutation_matrix(const struct inverse_matrix *m)
  91 {
  92         const unsigned *p = m->perm;
  93         const char *row[4] = {
  94                 "1 0 0 0\n",
  95                 "0 1 0 0\n",
  96                 "0 0 1 0\n",
  97                 "0 0 0 1\n"
  98         };
  99
 100         printf("  P =\n%s%s%s%s", row[p[0]], row[p[1]], row[p[2]], row[p[3]]);
 101 }
 102
 103 static void
 104 print_LU_decomposition(const struct inverse_matrix *m)
 105 {
 106         unsigned r, c;
 107
 108         printf("                            L                      "
 109                "                               U\n");
 110         for (r = 0; r < 4; ++r) {
 111                 double v;
 112
 113                 for (c = 0; c < 4; ++c) {
 114                         if (c < r)
 115                                 v = m->LU[r + c * 4];
 116                         else if (c == r)
 117                                 v = 1.0;
 118                         else
 119                                 v = 0.0;
 120                         printf(" %12.6f", v);
 121                 }
 122
 123                 printf(" | ");
 124
 125                 for (c = 0; c < 4; ++c) {
 126                         if (c >= r)
 127                                 v = m->LU[r + c * 4];
 128                         else
 129                                 v = 0.0;
 130                         printf(" %12.6f", v);
 131                 }
 132                 printf("\n");
 133         }
 134 }
 135
 136 static void
 137 print_inverse_data_matrix(const struct inverse_matrix *m)
 138 {
 139         unsigned r, c;
 140
 141         for (r = 0; r < 4; ++r) {
 142                 for (c = 0; c < 4; ++c)
 143                         printf(" %12.6f", m->LU[r + c * 4]);
 144                 printf("\n");
 145         }
 146
 147         printf("permutation: ");
 148         for (r = 0; r < 4; ++r)
 149                 printf(" %u", m->perm[r]);
 150         printf("\n");
 151 }
 152
 153 static void
 154 print_matrix(const struct weston_matrix *m)
 155 {
 156         unsigned r, c;
 157
 158         for (r = 0; r < 4; ++r) {
 159                 for (c = 0; c < 4; ++c)
 160                         printf(" %14.6e", m->d[r + c * 4]);
 161                 printf("\n");
 162         }
 163 }
 164
 165 static double
 166 frand(void)
 167 {
 168         double r = random();
 169         return r / (double)(RAND_MAX / 2) - 1.0f;
 170 }
 171
 172 static void
 173 randomize_matrix(struct weston_matrix *m)
 174 {
 175         unsigned i;
 176         for (i = 0; i < 16; ++i)
 177 #if 1
 178                 m->d[i] = frand() * exp(10.0 * frand());
 179 #else
 180                 m->d[i] = frand();
 181 #endif
 182 }
 183
 184 /* Take a matrix, compute inverse, multiply together
 185  * and subtract the identity matrix to get the error matrix.
 186  * Return the largest absolute value from the error matrix.
 187  */
 188 static double
 189 test_inverse(struct weston_matrix *m)
 190 {
 191         unsigned i;
 192         struct inverse_matrix q;
 193         double errsup = 0.0;
 194
 195         if (matrix_invert(q.LU, q.perm, m) != 0)
 196                 return INFINITY;
 197
 198         for (i = 0; i < 4; ++i)
 199                 inverse_transform(q.LU, q.perm, &m->d[i * 4]);
 200
 201         m->d[0] -= 1.0f;
 202         m->d[5] -= 1.0f;
 203         m->d[10] -= 1.0f;
 204         m->d[15] -= 1.0f;
 205
 206         for (i = 0; i < 16; ++i) {
 207                 double err = fabs(m->d[i]);
 208                 if (err > errsup)
 209                         errsup = err;
 210         }
 211
 212         return errsup;
 213 }
 214
 215 enum {
 216         TEST_OK,
 217         TEST_NOT_INVERTIBLE_OK,
 218         TEST_FAIL,
 219         TEST_COUNT
 220 };
 221
 222 static int
 223 test(void)
 224 {
 225         struct weston_matrix m;
 226         double det, errsup;
 227
 228         randomize_matrix(&m);
 229         det = determinant(&m);
 230
 231         errsup = test_inverse(&m);
 232         if (errsup < 1e-6)
 233                 return TEST_OK;
 234
 235         if (fabs(det) < 1e-5 && isinf(errsup))
 236                 return TEST_NOT_INVERTIBLE_OK;
 237
 238         printf("test fail, det: %g, error sup: %g\n", det, errsup);
 239
 240         return TEST_FAIL;
 241 }
 242
 243 static int running;
 244 static void
 245 stopme(int n)
 246 {
 247         running = 0;
 248 }
 249
 250 static void
 251 test_loop_precision(void)
 252 {
 253         int counts[TEST_COUNT] = { 0 };
 254
 255         printf("\nRunning a test loop for 10 seconds...\n");
 256         running = 1;
 257         alarm(10);
 258         while (running) {
 259                 counts[test()]++;
 260         }
 261
 262         printf("tests: %d ok, %d not invertible but ok, %d failed.\n"
 263                "Total: %d iterations.\n",
 264                counts[TEST_OK], counts[TEST_NOT_INVERTIBLE_OK],
 265                counts[TEST_FAIL],
 266                counts[TEST_OK] + counts[TEST_NOT_INVERTIBLE_OK] +
 267                counts[TEST_FAIL]);
 268 }
 269
 270 static void __attribute__((noinline))
 271 test_loop_speed_matrixvector(void)
 272 {
 273         struct weston_matrix m;
 274         struct weston_vector v = { { 0.5, 0.5, 0.5, 1.0 } };
 275         unsigned long count = 0;
 276         double t;
 277
 278         printf("\nRunning 3 s test on weston_matrix_transform()...\n");
 279
 280         weston_matrix_init(&m);
 281
 282         running = 1;
 283         alarm(3);
 284         reset_timer();
 285         while (running) {
 286                 weston_matrix_transform(&m, &v);
 287                 count++;
 288         }
 289         t = read_timer();
 290
 291         printf("%lu iterations in %f seconds, avg. %.1f us/iter.\n",
 292                count, t, 1e9 * t / count);
 293 }
 294
 295 static void __attribute__((noinline))
 296 test_loop_speed_inversetransform(void)
 297 {
 298         struct weston_matrix m;
 299         struct inverse_matrix inv;
 300         struct weston_vector v = { { 0.5, 0.5, 0.5, 1.0 } };
 301         unsigned long count = 0;
 302         double t;
 303
 304         printf("\nRunning 3 s test on inverse_transform()...\n");
 305
 306         weston_matrix_init(&m);
 307         matrix_invert(inv.LU, inv.perm, &m);
 308
 309         running = 1;
 310         alarm(3);
 311         reset_timer();
 312         while (running) {
 313                 inverse_transform(inv.LU, inv.perm, v.f);
 314                 count++;
 315         }
 316         t = read_timer();
 317
 318         printf("%lu iterations in %f seconds, avg. %.1f us/iter.\n",
 319                count, t, 1e9 * t / count);
 320 }
 321
 322 static void __attribute__((noinline))
 323 test_loop_speed_invert(void)
 324 {
 325         struct weston_matrix m;
 326         struct inverse_matrix inv;
 327         unsigned long count = 0;
 328         double t;
 329
 330         printf("\nRunning 3 s test on matrix_invert()...\n");
 331
 332         weston_matrix_init(&m);
 333
 334         running = 1;
 335         alarm(3);
 336         reset_timer();
 337         while (running) {
 338                 matrix_invert(inv.LU, inv.perm, &m);
 339                 count++;
 340         }
 341         t = read_timer();
 342
 343         printf("%lu iterations in %f seconds, avg. %.1f ns/iter.\n",
 344                count, t, 1e9 * t / count);
 345 }
 346
 347 static void __attribute__((noinline))
 348 test_loop_speed_invert_explicit(void)
 349 {
 350         struct weston_matrix m;
 351         unsigned long count = 0;
 352         double t;
 353
 354         printf("\nRunning 3 s test on weston_matrix_invert()...\n");
 355
 356         weston_matrix_init(&m);
 357
 358         running = 1;
 359         alarm(3);
 360         reset_timer();
 361         while (running) {
 362                 weston_matrix_invert(&m, &m);
 363                 count++;
 364         }
 365         t = read_timer();
 366
 367         printf("%lu iterations in %f seconds, avg. %.1f ns/iter.\n",
 368                count, t, 1e9 * t / count);
 369 }
 370
 371 int main(void)
 372 {
 373         struct sigaction ding;
 374         struct weston_matrix M;
 375         struct inverse_matrix Q;
 376         int ret;
 377         double errsup;
 378         double det;
 379
 380         ding.sa_handler = stopme;
 381         sigemptyset(&ding.sa_mask);
 382         ding.sa_flags = 0;
 383         sigaction(SIGALRM, &ding, NULL);
 384
 385         srandom(13);
 386
 387         M.d[0] = 3.0;   M.d[4] = 17.0;  M.d[8] = 10.0;  M.d[12] = 0.0;
 388         M.d[1] = 2.0;   M.d[5] = 4.0;   M.d[9] = -2.0;  M.d[13] = 0.0;
 389         M.d[2] = 6.0;   M.d[6] = 18.0;  M.d[10] = -12;  M.d[14] = 0.0;
 390         M.d[3] = 0.0;   M.d[7] = 0.0;   M.d[11] = 0.0;  M.d[15] = 1.0;
 391
 392         ret = matrix_invert(Q.LU, Q.perm, &M);
 393         printf("ret = %d\n", ret);
 394         printf("det = %g\n\n", determinant(&M));
 395
 396         if (ret != 0)
 397                 return 1;
 398
 399         print_inverse_data_matrix(&Q);
 400         printf("P * A = L * U\n");
 401         print_permutation_matrix(&Q);
 402         print_LU_decomposition(&Q);
 403
 404
 405         printf("a random matrix:\n");
 406         randomize_matrix(&M);
 407         det = determinant(&M);
 408         print_matrix(&M);
 409         errsup = test_inverse(&M);
 410         printf("\nThe matrix multiplied by its inverse, error:\n");
 411         print_matrix(&M);
 412         printf("max abs error: %g, original determinant %g\n", errsup, det);
 413
 414         test_loop_precision();
 415         test_loop_speed_matrixvector();
 416         test_loop_speed_inversetransform();
 417         test_loop_speed_invert();
 418         test_loop_speed_invert_explicit();
 419
 420         return 0;
 421 }