Devuan fork of gpsd
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

89 lines
2.3 KiB

  1. /*
  2. * Unit test for matrix-algebra code
  3. *
  4. * Check examples computed at
  5. * http://www.elektro-energetika.cz/calculations/matreg.php
  6. *
  7. * This file is Copyright (c) 2010 by the GPSD project
  8. * SPDX-License-Identifier: BSD-2-clause
  9. */
  10. #include <stdlib.h>
  11. #include <stdbool.h>
  12. #include <stdlib.h>
  13. #include <stdio.h>
  14. #include <math.h>
  15. #include "../compiler.h"
  16. #include "../matrix.h"
  17. static struct {
  18. double mat[4][4];
  19. double inv[4][4];
  20. } inverses[] = {
  21. /* identity matrix is self-inverse */
  22. {
  23. .mat = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}},
  24. .inv = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}},
  25. },
  26. /* inverse of a diagonal matrix has reciprocal values */
  27. {
  28. .mat = {{10,0,0,0}, {0,10,0,0}, {0,0,10,0}, {0,0,0,10}},
  29. .inv = {{0.1,0,0,0}, {0,0.1,0,0}, {0,0,0.1,0}, {0,0,0,0.1}},
  30. },
  31. /* random values with asymmetrical off-diagonal elements */
  32. {
  33. .mat = {{1,0,0,0}, {0,1,0,-2}, {0, 2,1,-4}, {0,0,0,1}},
  34. .inv = {{1,0,0,0}, {0,1,0, 2}, {0,-2,1, 0}, {0,0,0,1}},
  35. },
  36. {
  37. .mat = {{6,-4,1,-3},{-4,7,3,2},{1,3,6,-4},{-3,2,-4,6}},
  38. .inv = {{14,34,-40,-31},{34,84,-99,-77},{-40,-99,117,91},{-31,-77,91,71}},
  39. },
  40. };
  41. static void dump(const char *label, double m[4][4])
  42. {
  43. printf("%s:\n", label);
  44. printf("%f, %f, %f, %f\n", m[0][0], m[0][1], m[0][2], m[0][3]);
  45. printf("%f, %f, %f, %f\n", m[1][0], m[1][1], m[1][2], m[1][3]);
  46. printf("%f, %f, %f, %f\n", m[2][0], m[2][1], m[2][2], m[2][3]);
  47. printf("%f, %f, %f, %f\n", m[3][0], m[3][1], m[3][2], m[3][3]);
  48. }
  49. static bool check_diag(int n, double a[4][4], double b[4][4])
  50. {
  51. #define approx(x, y) (fabs(x - y) < 0.0001)
  52. if (approx(b[0][0], a[0][0]) && approx(b[1][1], a[1][1]) &&
  53. approx(b[2][2], a[2][2]) && approx(b[3][3], a[3][3]))
  54. return true;
  55. dump("a", a);
  56. dump("b", b);
  57. printf("Test %d residuals: %f %f %f %f\n",
  58. n,
  59. b[0][0] - a[0][0],
  60. b[1][1] - a[1][1],
  61. b[2][2] - a[2][2],
  62. b[3][3] - a[3][3]
  63. );
  64. return true;
  65. }
  66. int main(int argc UNUSED, char *argv[] UNUSED)
  67. {
  68. unsigned int i;
  69. for (i = 0; i < sizeof(inverses) / sizeof(inverses[0]); i++) {
  70. double inverse[4][4];
  71. if (!matrix_invert(inverses[i].mat, inverse)) {
  72. printf("Vanishing determinant in test %u\n", i);
  73. }
  74. if (!check_diag(i, inverse, inverses[i].inv))
  75. break;
  76. }
  77. printf("Matrix-algebra regression test succeeded\n");
  78. exit(0);
  79. }