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