1 //
2 // matrix_test.c
3 //
4 // This example tests basic matrix operations.
5 //
6
7 #include <stdio.h>
8 #include <string.h>
9 #include <math.h>
10 #include "liquid.h"
11
main()12 int main() {
13
14 float x[6] = {
15 1, 2, 3,
16 4, 5, 6};
17
18 float y[9] = {
19 1, 2, 3,
20 4, 5, 6,
21 7, 8, 9};
22
23 float z[6];
24
25 // compute z = x * y
26 printf("z = x * y :\n");
27 matrixf_mul(x,2,3,y,3,3,z,2,3);
28 matrixf_print(z,2,3);
29
30 /*
31 // compute z = y * x'
32 matrixf_transpose(x);
33 printf("x' : \n");
34 matrixf_print(x);
35 matrixf_transpose(z);
36 matrixf_multiply(y,x,z);
37 printf("z = y * x' :\n");
38 matrixf_print(z);
39
40 matrixf_destroy(x);
41 matrixf_destroy(y);
42 matrixf_destroy(z);
43 */
44
45 float s[16] = {
46 1,2,3,4,
47 5,5,7,8,
48 6,4,8,7,
49 1,0,3,1};
50 float s_inv[16];
51 memmove(s_inv,s,16*sizeof(float));
52 matrixf_inv(s_inv,4,4);
53
54 float i4[16];
55 matrixf_mul(s,4,4,s_inv,4,4,i4,4,4);
56
57 printf("\ns:\n");
58 matrixf_print(s,4,4);
59 printf("\ninv(s):\n");
60 matrixf_print(s_inv,4,4);
61 printf("\ns*inv(s):\n");
62 matrixf_print(i4,4,4);
63
64 printf("\n");
65 float det = matrixf_det(s,4,4);
66 printf("det(s) = %12.8f\n", det);
67
68 #if 0
69 // pivot test (matrix inversion)
70 float t[32] = {
71 1,2,3,4, 1,0,0,0,
72 5,5,7,8, 0,1,0,0,
73 6,4,8,7, 0,0,1,0,
74 1,0,3,1, 0,0,0,1};
75
76 unsigned int i;
77 for (i=0; i<4; i++) {
78 matrixf_pivot(t,4,8,i,i);
79 matrixf_print(t,4,8);
80 }
81
82 unsigned int j;
83 for (i=0; i<4; i++) {
84 float v = matrix_access(t,4,8,i,i);
85 for (j=0; j<8; j++)
86 matrix_access(t,4,8,i,j) /= v;
87 }
88 matrixf_print(t,4,8);
89 #endif
90
91 printf("\n");
92 printf("testing L/U decomposition [Crout's method]\n");
93 float L[16], U[16], P[16];
94 matrixf_ludecomp_crout(s,4,4,L,U,P);
95 matrixf_print(L,4,4);
96 matrixf_print(U,4,4);
97
98 printf("\n");
99 printf("testing L/U decomposition [Doolittle's method]\n");
100 matrixf_ludecomp_doolittle(s,4,4,L,U,P);
101 matrixf_print(L,4,4);
102 matrixf_print(U,4,4);
103
104 printf("\n\n");
105 float X[16] = {
106 0.84382, -2.38304, 1.43061, -1.66604,
107 3.99475, 0.88066, 4.69373, 0.44563,
108 7.28072, -2.06608, 0.67074, 9.80657,
109 6.07741, -3.93099, 1.22826, -0.42142};
110 float Y[16];
111 printf("\nX:\n");
112 matrixf_print(X,4,4);
113
114 // swaprows
115 memmove(Y,X,16*sizeof(float));
116 matrixf_swaprows(Y,4,4,0,2);
117 printf("\nmatrixf_swaprows(X,4,4,0,2):\n");
118 matrixf_print(Y,4,4);
119
120 // pivot test
121 memmove(Y,X,16*sizeof(float));
122 matrixf_pivot(Y,4,4,1,2);
123 printf("\nmatrixf_pivot(X,4,4,1,2):\n");
124 matrixf_print(Y,4,4);
125
126 // inverse test
127 memmove(Y,X,16*sizeof(float));
128 matrixf_inv(Y,4,4);
129 printf("\nmatrixf_inv(X,4,4):\n");
130 matrixf_print(Y,4,4);
131
132 // determinant test
133 float D = matrixf_det(X,4,4);
134 printf("\nmatrixf_det(X,4) = %12.8f\n", D);
135
136 // L/U decomp (Crout's method)
137 matrixf_ludecomp_crout(X,4,4,L,U,P);
138 printf("\nmatrixf_ludecomp_crout(X,4,4,L,U,P)\n");
139 printf("L:\n");
140 matrixf_print(L,4,4);
141 printf("U:\n");
142 matrixf_print(U,4,4);
143
144 // L/U decomp (Doolittle's method)
145 matrixf_ludecomp_doolittle(X,4,4,L,U,P);
146 printf("\nmatrixf_ludecomp_doolittle(X,4,4,L,U,P)\n");
147 printf("L:\n");
148 matrixf_print(L,4,4);
149 printf("U:\n");
150 matrixf_print(U,4,4);
151
152 printf("\n");
153 printf("testing Q/R decomposition [Gram-Schmidt]\n");
154 float Q[16], R[16];
155 matrixf_qrdecomp_gramschmidt(X,4,4,Q,R);
156 matrixf_print(Q,4,4);
157 matrixf_print(R,4,4);
158
159 /*
160 float b[4] = {
161 0.91489,
162 0.71789,
163 1.06553,
164 -0.81707};
165 */
166
167 float Xb[20] = {
168 0.84382, -2.38304, 1.43061, -1.66604, 0.91489,
169 3.99475, 0.88066, 4.69373, 0.44563, 0.71789,
170 7.28072, -2.06608, 0.67074, 9.80657, 1.06553,
171 6.07741, -3.93099, 1.22826, -0.42142, -0.81707};
172 printf("\n[X b] =\n");
173 matrixf_print(Xb,4,5);
174
175 matrixf_gjelim(Xb,4,5);
176 printf("\nmatrixf_gjelim(Xb,4,5)\n");
177 matrixf_print(Xb,4,5);
178
179 // compute a*a'
180 float a[20] = {
181 -0.24655, -1.78843, 0.39477, 0.43735, -1.08998,
182 -0.42751, 0.62496, 1.43802, 0.19814, 0.78155,
183 -0.35658, -0.81875, -1.09984, 1.87006, -0.94191,
184 0.39553, -2.02036, 1.17393, 1.54591, 1.29663};
185
186 printf("\na =\n");
187 matrixf_print(a,4,5);
188
189 printf("\n\n");
190 printf("computing a*a'\n");
191 float aaT[16];
192 matrixf_mul_transpose(a,4,5,aaT);
193 matrixf_print(aaT,4,4);
194
195 printf("\n\n");
196 printf("computing a'*a\n");
197 float aTa[25];
198 matrixf_transpose_mul(a,4,5,aTa);
199 matrixf_print(aTa,5,5);
200
201 printf("\n");
202 printf("testing Gram-Schmidt\n");
203 float Xgs[12] = {
204 1., 2., 1.,
205 0., 2., 0.,
206 2., 3., 1.,
207 1., 1., 0.};
208 float Ugs[12];
209 float Ugs_test[12] = {
210 sqrtf(6.)/6., sqrtf(2.)/6., 2./3.,
211 0., 2.*sqrtf(2.)/3.,-1./3.,
212 sqrtf(6.)/3., 0., 0.,
213 sqrtf(6.)/6., -sqrtf(2.)/6., -2./3.};
214 matrixf_gramschmidt(Xgs,4,3,Ugs);
215
216 printf("X:\n");
217 matrixf_print(Xgs,4,3);
218 printf("normalized X:\n");
219 matrixf_print(Ugs,4,3);
220 printf("expected:\n");
221 matrixf_print(Ugs_test,4,3);
222
223
224 //
225 // test Cholesky decomposition
226 //
227
228 printf("\n");
229 printf("testing Cholesky decomposition\n");
230 // generate input matrix
231 float complex Lp[9] = { 1.0, 0.0, 0.0,
232 -3.1 + 0.2*_Complex_I, 0.3, 0.0,
233 1.7 + 0.5*_Complex_I, -0.6 - 0.3*_Complex_I, 2.9};
234 float complex Ap[9];
235 matrixcf_mul_transpose(Lp, 3, 3, Ap);
236 float complex Lc[9];
237 matrixcf_chol(Ap, 3, Lc);
238
239 printf("Lp:\n"); matrixcf_print(Lp, 3, 3);
240 printf("Ap:\n"); matrixcf_print(Ap, 3, 3);
241 printf("Lc:\n"); matrixcf_print(Lc, 3, 3);
242
243 printf("done.\n");
244 return 0;
245 }
246
247