1 /*
2 * Copyright (c) 2019-2020 Kevin Townsend (KTOWN)
3 *
4 * SPDX-License-Identifier: Apache-2.0
5 */
6
7 #include <stdio.h>
8 #include <string.h>
9 #include <errno.h>
10 #include "zsl/matrices.h"
11 #include "zsl/vectors.h"
12
13 #define EPSILON 1e-4
14
svd_test()15 int svd_test()
16 {
17 size_t q = 18;
18 size_t p = 3;
19
20 zsl_real_t vi[3 * 18] = { 3.2, 5.0, -4.2, -3.1, 5.4, -2.2, 7.7, 8.7,
21 9.9, 8.9, 0.6, 5.4, 3.7, -3.3, 7.2, -5.4,
22 -0.6, 8.4, 2.4, 5.1, -5.7, 6.9, -2.1, 0.4,
23 6.4, -9.1, -4.4, 0.1, 7.4, 0.0, 6.1, -2.3,
24 5.5, 6.1, -8.4, -6.6, 7.1, 2.3, 4.1, -0.8,
25 -4.7, 9.9, -3.1, 1.2, 5.2, 6.4, -6.3, 5.2,
26 8.8, 7.3, 4.2, -4.7, 0.0, -0.4 };
27
28 /* Define the 18x3 matrix, assigning the values from 'vi'. */
29 ZSL_MATRIX_DEF(mep, q, p);
30 zsl_mtx_from_arr(&mep, vi);
31
32 /* Define the three SVD output matrices. */
33 ZSL_MATRIX_DEF(u, mep.sz_rows, mep.sz_rows);
34 ZSL_MATRIX_DEF(e, mep.sz_rows, mep.sz_cols);
35 ZSL_MATRIX_DEF(v, mep.sz_cols, mep.sz_cols);
36
37 /* Define placeholder matrices for SVD and related operations. */
38 ZSL_MATRIX_DEF(ue, mep.sz_rows, mep.sz_cols);
39 ZSL_MATRIX_DEF(vt, mep.sz_cols, mep.sz_cols);
40 ZSL_MATRIX_DEF(uevt, mep.sz_rows, mep.sz_cols);
41
42 /* Display the input matrix for reference. */
43 printf("INPUT (18x3 MATRIX)\n");
44 printf("-------------------\n");
45 zsl_mtx_print(&mep);
46 printf("\n");
47
48 /* Calculate and display the SVD using 150 iterations. */
49 zsl_mtx_svd(&mep, &u, &e, &v, 150);
50 printf("SVD OUTPUT\n");
51 printf("----------\n");
52 printf("U matrix (18x18):\n");
53 zsl_mtx_print(&u);
54 printf("\n");
55 printf("Sigma matrix (18x3):\n");
56 zsl_mtx_print(&e);
57 printf("\n");
58 printf("V matrix: (3x3)\n");
59 zsl_mtx_print(&v);
60 printf("\n");
61
62 /* Multiply u, sigma and v(transposed) from the SVD calculations,
63 * which shoud return the original matrix. */
64 printf("VERIFICATION\n");
65 printf("------------\n");
66 zsl_mtx_mult(&u, &e, &ue);
67 zsl_mtx_trans(&v, &vt);
68 zsl_mtx_mult(&ue, &vt, &uevt);
69 printf("U * Sigma * V^t:\n");
70 zsl_mtx_print(&uevt);
71 printf("\n");
72
73 return 0;
74 }
75
pinv_test()76 int pinv_test()
77 {
78 size_t q = 18;
79 size_t p = 3;
80
81 zsl_real_t vi[3 * 18] = { 3.2, 5.0, -4.2, -3.1, 5.4, -2.2, 7.7, 8.7,
82 9.9, 8.9, 0.6, 5.4, 3.7, -3.3, 7.2, -5.4,
83 -0.6, 8.4, 2.4, 5.1, -5.7, 6.9, -2.1, 0.4,
84 6.4, -9.1, -4.4, 0.1, 7.4, 0.0, 6.1, -2.3,
85 5.5, 6.1, -8.4, -6.6, 7.1, 2.3, 4.1, -0.8,
86 -4.7, 9.9, -3.1, 1.2, 5.2, 6.4, -6.3, 5.2,
87 8.8, 7.3, 4.2, -4.7, 0.0, -0.4 };
88
89 /* Define the 18x3 matrix, assigning the chosen values. */
90 ZSL_MATRIX_DEF(mep, q, p);
91 zsl_mtx_from_arr(&mep, vi);
92
93 /* Define the pseudoinverse output matrix. */
94 ZSL_MATRIX_DEF(pinv, p, q);
95
96 /* Display the input matrix for reference. */
97 printf("INPUT (18x3 MATRIX)\n");
98 printf("-------------------\n");
99 zsl_mtx_print(&mep);
100 printf("\n");
101
102 /* Calculate and display the pseudoinverse using 150 iterations. */
103 zsl_mtx_pinv(&mep, &pinv, 150);
104 printf("PSEUDOINVERSE (3x18 MATRIX)\n");
105 printf("--------------------------\n");
106 zsl_mtx_print(&pinv);
107 printf("\n");
108
109 return 0;
110 }
111
112 int
main(void)113 main(void)
114 {
115 printf("Hello, zscilib!\n\n");
116
117 printf("SVD TEST\n");
118 printf("--------\n");
119 svd_test();
120
121 printf("PSEUDOINVERSE TEST\n");
122 printf("------------------\n");
123 pinv_test();
124
125 return 0;
126 }
127
