check_linear_algebra.c
Go to the documentation of this file.
00001 #include <math.h>
00002 #include <stdlib.h>
00003 #include <check.h>
00004 
00005 #include <stdio.h>
00006 #include "check_utils.h"
00007 
00008 #include <linear_algebra.h>
00009 
00010 #define LINALG_TOL 1e-10
00011 #define LINALG_NUM 222
00012 #define MATRIX_MIN -1e22
00013 #define MATRIX_MAX 1e22
00014 #define mrand frand(MATRIX_MIN, MATRIX_MAX)
00015 #define MSIZE_MAX 64
00016 
00017 #define VEC_PRINTF(v, _n) {                                     \
00018     printf("%s:%u <|%s| %lf",                                   \
00019            __FILE__, __LINE__, #v, v[0]);                       \
00020     for (int _i = 1; _i < _n; _i++) printf(", %lf", v[_i]);     \
00021     printf(">\n");                                              \
00022   }
00023 
00024 #define MAT_PRINTF(m, _r, _c) {                    \
00025     printf("%s:%u <|%s|",                          \
00026            __FILE__, __LINE__, #m);                \
00027     for (int _i = 0; _i < _r; _i++) {              \
00028       printf(" [%lf", m[_i*_c + 0]);                \
00029       for (int _j = 1; _j < _c; _j++)              \
00030         printf(" %lf", m[_i*_c + _j]);               \
00031       printf("]");                                 \
00032       if (_r > 2 && _i < _r - 1) printf("\n\t\t\t       ");              \
00033     }                                              \
00034     printf(">\n");                                 \
00035   }
00036 
00037 /* TODO: matrix_multiply, matrix_add_sc, matrix_copy, all vector functions */
00038 
00039 START_TEST(test_matrix_inverse_2x2) {
00040   u32 i, j, t;
00041   double A[4];
00042   double B[4];
00043   double I[4];
00044 
00045   seed_rng();
00046   /* 2x2 inverses */
00047   for (t = 0; t < LINALG_NUM; t++) {
00048     do {
00049       for (i = 0; i < 2; i++)
00050         for (j = 0; j < 2; j++)
00051           A[2*i + j] = mrand;
00052     } while (matrix_inverse(2, A, B) < 0);
00053     matrix_multiply(2, 2, 2, A, B, I);
00054     fail_unless(fabs(I[0] - 1) < LINALG_TOL,
00055                 "Matrix differs from identity: %lf", I[0]);
00056     fail_unless(fabs(I[3] - 1) < LINALG_TOL,
00057                 "Matrix differs from identity: %lf", I[3]);
00058   }
00059   for (i = 0; i < 2; i++)
00060     for (j = 0; j < 2; j++)
00061       if (j == 0)
00062         A[2*i + j] = 22;
00063       else
00064         A[2*i + j] = 1;
00065   s32 mi = matrix_inverse(2, A, B);
00066   fail_unless(mi < 0, "Singular matrix not detected.");
00067 }
00068 END_TEST
00069 
00070 START_TEST(test_matrix_inverse_3x3) {
00071   u32 i, j, t;
00072   double A[9];
00073   double B[9];
00074   double I[9];
00075 
00076   seed_rng();
00077   /* 3x3 inverses */
00078   for (t = 0; t < LINALG_NUM; t++) {
00079     do {
00080       for (i = 0; i < 3; i++)
00081         for (j = 0; j < 3; j++)
00082           A[3*i + j] = mrand;
00083     } while (matrix_inverse(3, A, B) < 0);
00084     matrix_multiply(3, 3, 3, A, B, I);
00085     fail_unless(fabs(I[0] - 1) < LINALG_TOL,
00086                 "Matrix differs from identity: %lf", I[0]);
00087     fail_unless(fabs(I[4] - 1) < LINALG_TOL,
00088                 "Matrix differs from identity: %lf", I[4]);
00089     fail_unless(fabs(I[8] - 1) < LINALG_TOL,
00090                 "Matrix differs from identity: %lf", I[8]);
00091   }
00092   for (i = 0; i < 3; i++)
00093     for (j = 0; j < 3; j++)
00094       if (j == 0)
00095         A[3*i + j] = 33;
00096       else
00097         A[3*i + j] = 1;
00098   s32 mi = matrix_inverse(3, A, B);
00099   fail_unless(mi < 0, "Singular matrix not detected.");
00100 }
00101 END_TEST
00102 
00103 START_TEST(test_matrix_inverse_4x4) {
00104   u32 i, j, t;
00105   double A[16];
00106   double B[16];
00107   double I[16];
00108 
00109   seed_rng();
00110   /* 4x4 inverses */
00111   for (t = 0; t < LINALG_NUM; t++) {
00112     do {
00113       for (i = 0; i < 4; i++)
00114         for (j = 0; j < 4; j++)
00115           A[4*i + j] = mrand;
00116     } while (matrix_inverse(4, A, B) < 0);
00117     matrix_multiply(4, 4, 4, A, B, I);
00118     fail_unless(fabs(I[0] - 1) < LINALG_TOL,
00119                 "Matrix differs from identity: %lf", I[0]);
00120     fail_unless(fabs(I[5] - 1) < LINALG_TOL,
00121                 "Matrix differs from identity: %lf", I[5]);
00122     fail_unless(fabs(I[10] - 1) < LINALG_TOL,
00123                 "Matrix differs from identity: %lf", I[10]);
00124     fail_unless(fabs(I[15] - 1) < LINALG_TOL,
00125                 "Matrix differs from identity: %lf", I[15]);
00126   }
00127   for (i = 0; i < 4; i++)
00128     for (j = 0; j < 4; j++)
00129       if (j == 0)
00130         A[4*i + j] = 44;
00131       else
00132         A[4*i + j] = 1;
00133   s32 mi = matrix_inverse(4, A, B);
00134   fail_unless(mi < 0, "Singular matrix not detected.");
00135 }
00136 END_TEST
00137 
00138 START_TEST(test_matrix_inverse_5x5) {
00139   u32 i, j, t;
00140   double A[25];
00141   double B[25];
00142   double I[25];
00143 
00144   seed_rng();
00145   /* 5x5 inverses */
00146   for (t = 0; t < LINALG_NUM; t++) {
00147     do {
00148       for (i = 0; i < 5; i++)
00149         for (j = 0; j < 5; j++)
00150           A[5*i + j] = mrand;
00151     } while (matrix_inverse(5, A, B) < 0);
00152     matrix_multiply(5, 5, 5, A, B, I);
00153     fail_unless(fabs(I[0] - 1) < LINALG_TOL,
00154                 "Matrix differs from identity: %lf", I[0]);
00155     fail_unless(fabs(I[6] - 1) < LINALG_TOL,
00156                 "Matrix differs from identity: %lf", I[6]);
00157     fail_unless(fabs(I[12] - 1) < LINALG_TOL,
00158                 "Matrix differs from identity: %lf", I[12]);
00159     fail_unless(fabs(I[18] - 1) < LINALG_TOL,
00160                 "Matrix differs from identity: %lf", I[18]);
00161     fail_unless(fabs(I[24] - 1) < LINALG_TOL,
00162                 "Matrix differs from identity: %lf", I[24]);
00163   }
00164   for (i = 0; i < 5; i++)
00165     for (j = 0; j < 5; j++)
00166       if (j == 0)
00167         A[5*i + j] = 55;
00168       else
00169         A[5*i + j] = 1;
00170   s32 mi = matrix_inverse(5, A, B);
00171   fail_unless(mi < 0, "Singular matrix not detected.");
00172 }
00173 END_TEST
00174 
00175 START_TEST(test_matrix_add_sc) {
00176   u32 i, j, t;
00177 
00178   seed_rng();
00179   for (t = 0; t < LINALG_NUM; t++) {
00180     u32 n = sizerand(MSIZE_MAX);
00181     u32 m = sizerand(MSIZE_MAX);
00182     double A[n*m];
00183     double B[m*n];
00184     for (i = 0; i < n; i++)
00185       for (j = 0; j < m; j++) {
00186         A[m*i + j] = mrand;
00187       }
00188     matrix_add_sc(n, m, A, A, -1, B);
00189     for (i = 0; i < n; i++)
00190       for (j = 0; j < m; j++)
00191         fail_unless(fabs(B[m*i + j]) < LINALG_TOL,
00192                     "Matrix differs from zero: %lf", B[m*i + j]);
00193   }
00194 }
00195 END_TEST
00196 
00197 START_TEST(test_matrix_copy) {
00198   u32 i, j, t;
00199   double tmp;
00200 
00201   seed_rng();
00202   for (t = 0; t < LINALG_NUM; t++) {
00203     u32 n = sizerand(MSIZE_MAX);
00204     u32 m = sizerand(MSIZE_MAX);
00205     double A[n*m];
00206     double B[m*n];
00207     for (i = 0; i < n; i++)
00208       for (j = 0; j < m; j++)
00209         A[m*i + j] = mrand;
00210     matrix_copy(n, m, A, B);
00211     for (i = 0; i < n; i++)
00212       for (j = 0; j < m; j++) {
00213         tmp = fabs(B[m*i + j] - A[m*i + j]);
00214         fail_unless(tmp < LINALG_TOL,
00215                     "Matrix differs from zero: %lf", tmp);
00216       }
00217   }
00218 }
00219 END_TEST
00220 
00221 START_TEST(test_matrix_transpose) {
00222   u32 i, j, t;
00223 
00224   seed_rng();
00225   for (t = 0; t < LINALG_NUM; t++) {
00226     u32 n = sizerand(MSIZE_MAX);
00227     u32 m = sizerand(MSIZE_MAX);
00228     double A[n*m];
00229     double B[m*n];
00230     double C[n*m];
00231     for (i = 0; i < n; i++)
00232       for (j = 0; j < m; j++)
00233         A[m*i + j] = mrand;
00234     matrix_transpose(n, m, A, B);
00235     matrix_transpose(m, n, B, C);
00236     for (i = 0; i < n; i++)
00237       for (j = 0; j < m; j++)
00238         fail_unless(fabs(A[m*i + j] - C[m*i + j]) < LINALG_TOL,
00239                     "Matrix element differs from original: %lf, %lf",
00240                     A[m*i + j], C[m*i + j]);
00241   }
00242 }
00243 END_TEST
00244 
00245 START_TEST(test_vector_dot) {
00246   u32 i, t;
00247 
00248   seed_rng();
00249   for (t = 0; t < LINALG_NUM; t++) {
00250     u32 n = sizerand(MSIZE_MAX);
00251     u32 mid;
00252     if (n % 2 == 0)
00253       mid = n / 2;
00254     else
00255       mid = (n - 1) / 2;
00256 
00257     double A[n], B[n];
00258     for (i = 0; i < n; i++) {
00259       A[i] = mrand/1e20;
00260       if (i < mid)
00261         B[n-i-1] = -A[i];
00262       else
00263         B[n-i-1] = A[i];
00264     }
00265     double dot = vector_dot(n, A, B);
00266     if (n % 2 == 0)
00267       fail_unless(fabs(dot) < LINALG_TOL,
00268                   "Dot product differs from zero: %lf", vector_dot(n, A, B));
00269     else
00270       fail_unless(fabs(dot - A[mid]*B[mid]) < LINALG_TOL,
00271                   "Dot product differs from square of middle element "
00272                   "%lf: %lf (%lf)",
00273                   A[mid]*B[mid], dot, dot - A[mid]*B[mid]);
00274   }
00275 }
00276 END_TEST
00277 
00278 START_TEST(test_vector_mean) {
00279   u32 i, t;
00280 
00281   seed_rng();
00282   for (t = 0; t < LINALG_NUM; t++) {
00283     u32 n = sizerand(MSIZE_MAX);
00284     double A[n];
00285     double test = mrand/1e22;
00286     for (i = 0; i < n; i++)
00287       A[i] = test + i;
00288     double mean = vector_mean(n, A);
00289     double expect = test + (n - 1.0) / 2.0;
00290     fail_unless(fabs(mean - expect)
00291                 < LINALG_TOL,
00292                 "Mean differs from expected %lf: %lf (%lf)",
00293                 expect, mean, fabs(mean - expect));
00294   }
00295 }
00296 END_TEST
00297 
00298 START_TEST(test_vector_norm) {
00299   u32 i, t;
00300 
00301   seed_rng();
00302   for (t = 0; t < LINALG_NUM; t++) {
00303     u32 n = sizerand(MSIZE_MAX);
00304     double test = mrand/1e22;
00305     double A[n];
00306     for (i = 0; i < n; i++)
00307       A[i] = test;
00308     fail_unless(fabs(vector_norm(n, A)*vector_norm(n, A) - n * test * test)
00309                 < LINALG_TOL * vector_norm(n, A),
00310                 "Norm differs from expected %lf: %lf (%lf)",
00311                 n * test * test, vector_norm(n, A) * vector_norm(n, A),
00312                 fabs(vector_norm(n, A) * vector_norm(n, A) - n * test * test));
00313   }
00314 }
00315 END_TEST
00316 
00317 START_TEST(test_vector_normalize) {
00318   u32 i, t;
00319 
00320   seed_rng();
00321   for (t = 0; t < LINALG_NUM; t++) {
00322     u32 n = sizerand(MSIZE_MAX);
00323     double A[n];
00324     for (i = 0; i < n; i++)
00325       A[i] = mrand;
00326     vector_normalize(n, A);
00327     double vnorm = vector_norm(n, A);
00328     fail_unless(fabs(vnorm - 1) < LINALG_TOL,
00329                 "Norm differs from 1: %lf",
00330                 vector_norm(n, A));
00331   }
00332 }
00333 END_TEST
00334 
00335 START_TEST(test_vector_add_sc) {
00336   u32 i, t;
00337 
00338   seed_rng();
00339   for (t = 0; t < LINALG_NUM; t++) {
00340     u32 n = sizerand(MSIZE_MAX);
00341     double A[n], B[n];
00342     for (i = 0; i < n; i++)
00343       A[i] = mrand;
00344     vector_add_sc(n, A, A, -1, B);
00345     for (i = 0; i < n; i++)
00346       fail_unless(fabs(B[i]) < LINALG_TOL,
00347                   "Vector element differs from 0: %lf",
00348                   B[i]);
00349   }
00350 }
00351 END_TEST
00352 
00353 START_TEST(test_vector_add) {
00354   u32 i, t;
00355 
00356   seed_rng();
00357   for (t = 0; t < LINALG_NUM; t++) {
00358     u32 n = sizerand(MSIZE_MAX);
00359     double A[n], B[n], C[n];
00360     for (i = 0; i < n; i++) {
00361       A[i] = mrand;
00362       B[i] = -A[i];
00363     }
00364     vector_add(n, A, B, C);
00365     for (i = 0; i < n; i++)
00366       fail_unless(fabs(C[i]) < LINALG_TOL,
00367                   "Vector element differs from 0: %lf",
00368                   C[i]);
00369   }
00370 }
00371 END_TEST
00372 
00373 START_TEST(test_vector_subtract) {
00374   u32 i, t;
00375 
00376   seed_rng();
00377   for (t = 0; t < LINALG_NUM; t++) {
00378     u32 n = sizerand(MSIZE_MAX);
00379     double A[n], B[n], C[n];
00380     for (i = 0; i < n; i++) {
00381       A[i] = mrand;
00382       B[i] = A[i];
00383     }
00384     vector_subtract(n, A, B, C);
00385     for (i = 0; i < n; i++)
00386       fail_unless(fabs(C[i]) < LINALG_TOL,
00387                   "Vector element differs from 0: %lf",
00388                   C[i]);
00389   }
00390 }
00391 END_TEST
00392 
00393 START_TEST(test_vector_cross) {
00394   u32 i, t;
00395 
00396   seed_rng();
00397   for (t = 0; t < LINALG_NUM; t++) {
00398     double A[3], B[3], C[3], D[3];
00399     for (i = 0; i < 3; i++) {
00400       A[i] = mrand;
00401       B[i] = A[i];
00402     }
00403     vector_cross(A, B, C);
00404     for (i = 0; i < 3; i++)
00405       fail_unless(fabs(C[i]) < LINALG_TOL,
00406                   "Vector element differs from 0: %lf",
00407                   C[i]);
00408     for (i = 0; i < 3; i++) {
00409       A[i] = mrand;
00410       B[i] = mrand;
00411     }
00412     vector_cross(A, B, C);
00413     for (i = 0; i < 3; i++) B[i] *= -1;
00414     vector_cross(B, A, D);
00415     for (i = 0; i < 3; i++)
00416       fail_unless(fabs(C[i] - D[i]) < LINALG_TOL,
00417                   "Vector equality fails: %lf != %lf",
00418                   C[i], D[i]);
00419   }
00420 }
00421 END_TEST
00422 
00423 START_TEST(test_vector_three) {
00424   u32 i, t;
00425 
00426   seed_rng();
00427   for (t = 0; t < LINALG_NUM; t++) {
00428     double A[3], B[3], C[3], tmp[3];
00429     double D, E, F, norm;
00430     for (i = 0; i < 3; i++) {
00431       A[i] = mrand / 1e20;
00432       B[i] = mrand / 1e20;
00433       C[i] = mrand / 1e20;
00434     }
00435     /* Check triple product identity */
00436     vector_cross(A, B, tmp);
00437     D = vector_dot(3, C, tmp);
00438     vector_cross(B, C, tmp);
00439     E = vector_dot(3, A, tmp);
00440     vector_cross(C, A, tmp);
00441     F = vector_dot(3, B, tmp);
00442 
00443     norm = (vector_norm(3, A) + vector_norm(3, B) + vector_norm(3, C))/3;
00444     fail_unless(fabs(E - D) < LINALG_TOL * norm,
00445                 "Triple product failure between %lf and %lf",
00446                 D, E);
00447     fail_unless(fabs(E - F) < LINALG_TOL * norm,
00448                 "Triple product failure between %lf and %lf",
00449                 E, F);
00450     fail_unless(fabs(F - D) < LINALG_TOL * norm,
00451                 "Triple product failure between %lf and %lf",
00452                 F, D);
00453   }
00454 }
00455 END_TEST
00456 
00457 /*
00458 START_TEST(test_qrsolve_consistency) {
00459   u32 i, j, t;
00460   double norm;
00461 
00462   seed_rng();
00463   for (t = 0; t < LINALG_NUM; t++) {
00464     u32 n = sizerand(MSIZE_MAX);
00465     double x_gauss[n], x_qr[n];
00466     double A[n*n], Ainv[n*n], b[n];
00467     do {
00468       for (i = 0; i < n; i++) {
00469         b[i] = mrand;
00470         for (j = 0; j < n; j++)
00471           A[n*i + j] = mrand;
00472       }
00473     } while (matrix_inverse(n, A, Ainv) < 0);
00474     matrix_multiply(n, n, 1, Ainv, b, x_gauss);
00475     qrsolve(A, n, n, b, x_qr);
00476 
00477     norm = (vector_norm(n, x_qr) + vector_norm(n, x_gauss)) / 2;
00478     for (i = 0; i < n; i++)
00479       fail_unless(fabs(x_qr[i] - x_gauss[i]) < LINALG_TOL * norm,
00480                   "QR solve failure; difference was %lf for element %u",
00481                   x_qr[i] - x_gauss[i], i);
00482   }
00483 }
00484 END_TEST
00485 
00486 START_TEST(test_qrsolve_rect) {
00487   s32 i;
00488   const double A[8] = {-0.0178505395610981, 1.4638781031761146,
00489                        -0.8242742209580581, -0.6843477128009663,
00490                        0.9155272861151404, -0.1651159277864960,
00491                        -0.9929037180867774, -0.1491537478964264};
00492   double Q[16], R[8];
00493 
00494   double buf[10] __attribute__((unused)) = {22, 22, 22, 22, 22,
00495                                             22, 22, 22, 22, 22};
00496 
00497   i = qrdecomp(A, 4, 2, Q, R);
00498 
00499   printf("i returned %d\n", i);
00500 
00501   MAT_PRINTF(A, 4, 2);
00502   MAT_PRINTF(Q, 4, 4);
00503   MAT_PRINTF(R, 4, 2);
00504 }
00505 END_TEST
00506 */
00507 
00508 Suite* linear_algebra_suite(void) {
00509   Suite *s = suite_create("Linear algebra");
00510 
00511   /* Core test case */
00512   TCase *tc_core = tcase_create("Core");
00513   tcase_add_test(tc_core, test_matrix_add_sc);
00514   tcase_add_test(tc_core, test_matrix_copy);
00515   tcase_add_test(tc_core, test_matrix_transpose);
00516   tcase_add_test(tc_core, test_matrix_inverse_2x2);
00517   tcase_add_test(tc_core, test_matrix_inverse_3x3);
00518   tcase_add_test(tc_core, test_matrix_inverse_4x4);
00519   tcase_add_test(tc_core, test_matrix_inverse_5x5);
00520 
00521   tcase_add_test(tc_core, test_vector_dot);
00522   tcase_add_test(tc_core, test_vector_mean);
00523   tcase_add_test(tc_core, test_vector_norm);
00524   tcase_add_test(tc_core, test_vector_normalize);
00525   tcase_add_test(tc_core, test_vector_add_sc);
00526   tcase_add_test(tc_core, test_vector_add);
00527   tcase_add_test(tc_core, test_vector_subtract);
00528   tcase_add_test(tc_core, test_vector_cross);
00529   tcase_add_test(tc_core, test_vector_three);
00530   /*tcase_add_test(tc_core, test_qrsolve_consistency);*/
00531   /*tcase_add_test(tc_core, test_qrsolve_rect);*/
00532   suite_add_tcase(s, tc_core);
00533 
00534   return s;
00535 }
00536 


enu
Author(s): Mike Purvis
autogenerated on Fri Jan 3 2014 11:21:07