|
| 1 | +#include <stdio.h> |
| 2 | +#include <stdlib.h> |
| 3 | + |
| 4 | +#include <check.h> |
| 5 | + |
| 6 | +#include <libswiftnav/linear_algebra.h> |
| 7 | +#include <libswiftnav/printing_utils.h> |
| 8 | +#include <libswiftnav/plover/qr.h> |
| 9 | + |
| 10 | +#include "check_utils.h" |
| 11 | + |
| 12 | +START_TEST(test_ok) |
| 13 | +{ |
| 14 | + double m1[] = {1, 0, 2, |
| 15 | + 1, 1, -1, |
| 16 | + 0, 0, 1, |
| 17 | + 0, 0, 22 }; |
| 18 | + double b[] = {5, 4, 1.1, 22}; |
| 19 | + double soln[3]; |
| 20 | + double residual; |
| 21 | + |
| 22 | + s8 code = qr_solve(4, 3, m1, b, soln, &residual); |
| 23 | + |
| 24 | + fail_unless(code == 0, "Solver error code: %d\n", code); |
| 25 | + fail_unless(residual < 0.1, "Residual too large: %f\n", residual); |
| 26 | +} |
| 27 | +END_TEST |
| 28 | + |
| 29 | +/* Deficient */ |
| 30 | +START_TEST(test_bad) |
| 31 | +{ |
| 32 | + double m1[] = {1, 0, 0, |
| 33 | + 0, 1, 0, |
| 34 | + 0, 0, 0 }; |
| 35 | + double b[] = {1,1,1}; |
| 36 | + double soln[3]; |
| 37 | + double residual; |
| 38 | + |
| 39 | + s8 code = qr_solve(3, 3, m1, b, soln, &residual); |
| 40 | + |
| 41 | + fail_unless(code == -1, "Wrong solver code: %d\n", code); |
| 42 | +} |
| 43 | +END_TEST |
| 44 | + |
| 45 | +double randd() { |
| 46 | + int max = 22; |
| 47 | + return (double)(rand()) / RAND_MAX * max - (double)max / 2; |
| 48 | +} |
| 49 | + |
| 50 | +void test_qr(int rows, int cols) |
| 51 | +{ |
| 52 | + double mat[rows*cols]; |
| 53 | + double vec[cols]; |
| 54 | + |
| 55 | + for(int r = 0; r < rows; r++) { |
| 56 | + for(int c = 0; c < cols; c++) { |
| 57 | + mat[r*cols + c] = randd(); |
| 58 | + } |
| 59 | + } |
| 60 | + for(int c = 0; c < cols; c++) { |
| 61 | + vec[c] = randd(); |
| 62 | + } |
| 63 | + |
| 64 | + double m1[rows*cols]; |
| 65 | + memcpy(m1, mat, sizeof(m1)); |
| 66 | + |
| 67 | + //printf("%d %d\n", rows, cols); |
| 68 | + //print_double_mtx(mat, rows, cols); |
| 69 | + |
| 70 | + double out[rows]; |
| 71 | + matrix_multiply(rows, cols, 1, mat, vec, out); |
| 72 | + double soln[cols]; |
| 73 | + double residual; |
| 74 | + int code = qr_solve(rows, cols, mat, out, soln, &residual); |
| 75 | + fail_unless(residual < 0.0000001, "Residual too large: %f\n", residual); |
| 76 | + if (!arr_within_epsilon(cols, soln, vec) && code != -1) { |
| 77 | + printf("original matrix: "); |
| 78 | + print_double_mtx(m1, rows, cols); |
| 79 | + printf("reduced R matrix: "); |
| 80 | + print_double_mtx(mat, rows, cols); |
| 81 | + printf("solution: "); |
| 82 | + print_double_mtx(soln, 1, cols); |
| 83 | + printf("vec: "); |
| 84 | + print_double_mtx(vec, 1, cols); |
| 85 | + fail_unless(false, "Solution wrong\n"); |
| 86 | + } |
| 87 | +} |
| 88 | + |
| 89 | +/* Performs 10000 trials. */ |
| 90 | +START_TEST(test_rand) |
| 91 | +{ |
| 92 | + srand(0); |
| 93 | + for (int i = 0; i < 10000; i++) { |
| 94 | + int rows = rand() % 25 + 1; |
| 95 | + int cols = rand() % rows + 1; |
| 96 | + test_qr(rows, cols); |
| 97 | + } |
| 98 | +} |
| 99 | +END_TEST |
| 100 | + |
| 101 | +Suite* qr_test_suite(void) |
| 102 | +{ |
| 103 | + Suite *s = suite_create("Generated QR solver"); |
| 104 | + |
| 105 | + TCase *tc_core = tcase_create("Core"); |
| 106 | + tcase_add_test(tc_core, test_ok); |
| 107 | + tcase_add_test(tc_core, test_bad); |
| 108 | + tcase_add_test(tc_core, test_rand); |
| 109 | + suite_add_tcase(s, tc_core); |
| 110 | + return s; |
| 111 | +} |
0 commit comments