RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
ieee_mobile_computing_lin_alg.c
1 /*
2  * ieee_mobile_computing.c
3  *
4  * Created on: 26.09.2019
5  * Author: zakaria
6  */
7 
8 #include <stdio.h>
9 #include "matrix.h"
10 #include "vector.h"
11 #include "qr_givens.h"
12 #include "qr_householder.h"
13 #include "qr_common.h"
14 #include "svd.h"
15 #include "solve.h"
16 
17 #define SEC 1000000
18 #define ROW_MAX 10
19 #define COL_MAX 5
20 
21 matrix_t A[10][5] = { { 0.8147, 0.1576, 0.6557, 0.7060, 0.4387 },
22  { 0.9058, 0.9706, 0.0357, 0.0318, 0.3816 },
23  { 0.1270, 0.9572, 0.8491, 0.2769, 0.7655 },
24  { 0.9134, 0.4854, 0.9340, 0.0462, 0.7952 },
25  { 0.6324, 0.8003, 0.6787, 0.0971, 0.1869 },
26  { 0.0975, 0.1419, 0.7577, 0.8235, 0.4898 },
27  { 0.2785, 0.4218, 0.7431, 0.6948, 0.4456 },
28  { 0.5469, 0.9157, 0.3922, 0.3171, 0.6463 },
29  { 0.9575, 0.7922, 0.6555, 0.9502, 0.7094 },
30  { 0.9649, 0.9595, 0.1712, 0.0344, 0.7547 } };
31 
32 matrix_t copy_A[10][5] = { { 0.8147, 0.1576, 0.6557, 0.7060, 0.4387 },
33  { 0.9058, 0.9706, 0.0357, 0.0318, 0.3816 },
34  { 0.1270, 0.9572, 0.8491, 0.2769, 0.7655 },
35  { 0.9134, 0.4854, 0.9340, 0.0462, 0.7952 },
36  { 0.6324, 0.8003, 0.6787, 0.0971, 0.1869 },
37  { 0.0975, 0.1419, 0.7577, 0.8235, 0.4898 },
38  { 0.2785, 0.4218, 0.7431, 0.6948, 0.4456 },
39  { 0.5469, 0.9157, 0.3922, 0.3171, 0.6463 },
40  { 0.9575, 0.7922, 0.6555, 0.9502, 0.7094 },
41  { 0.9649, 0.9595, 0.1712, 0.0344, 0.7547 } };
42 
43 matrix_t B[5][10] = { { 0.2760, 0.4984, 0.7513, 0.9593, 0.8407, 0.3500, 0.3517, 0.2858, 0.0759, 0.1299 },
44  { 0.6797, 0.9597, 0.2551, 0.5472, 0.2543, 0.1966, 0.8308, 0.7572, 0.0540, 0.5688 },
45  { 0.6551, 0.3404, 0.5060, 0.1386, 0.8143, 0.2511, 0.5853, 0.7537, 0.5308, 0.4694 },
46  { 0.1626, 0.5853, 0.6991, 0.1493, 0.2435, 0.6160, 0.5497, 0.3804, 0.7792, 0.0119 },
47  { 0.1190, 0.2238, 0.8909, 0.2575, 0.9293, 0.4733, 0.9172, 0.5678, 0.9340, 0.3371 } };
48 
49 matrix_t b[10] = { 1.2902, 0.8819, 0.9721, 1.2347, 0.9185, 0.9844, 1.0627, 1.0280, 1.7283, 1.0618 };
50 
51 void ieee_mobile_comp_test_matrix_mul(void)
52 {
53  puts("########### Measures of the times of matrices multiplications ###########");
54  matrix_t C[10][10];
55 
56  for (uint8_t i = 0; i < ROW_MAX; i++) {
57  // Start timer
58  matrix_part_mul(COL_MAX, A, ROW_MAX, B, 0, i, 0, COL_MAX - 1, 0, COL_MAX - 1, 0, ROW_MAX - 1, ROW_MAX, C);
59  // Stop timer
60 
61  printf("C_%u = ", i + 1);
62  matrix_flex_part_print(ROW_MAX, ROW_MAX, C, 0, i, 0, ROW_MAX - 1, 1, 4);
63  puts(" ");
64  }
65 }
66 
67 void ieee_solve_lin_equ_sys(void)
68 {
69 
70  puts("Square sub-matrices:");
71 
72  for (uint8_t i = 1; i < COL_MAX; i++) {
73  matrix_t sub_A[i + 1][i + 1];
74  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
75  i + 1, sub_A);
76  printf("(%u,%u) matrix: ", i + 1, i + 1);
77  matrix_flex_print(i + 1, i + 1, sub_A, 3, 4);
78  puts("");
79  enum ALGORITHM algo = Moore_Penrose;
80  puts("Moore Penrose: ");
81  matrix_t x[i + 1];
82 
83  vector_clear(i+1, x);
84  // Start Timer !!!
85  int8_t error_no = solve(i + 1, i + 1, sub_A, b, x, algo);
86  // Stop Timer !!!
87 
88  if (error_no > 0) {
89  printf("x = ");
90  vector_flex_print(i + 1, x, 3, 4);
91  puts("");
92  }
93  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
94  i + 1, sub_A);
95  vector_t A_mul_x[i+1];
96  matrix_mul_vec(i+1, i+1, sub_A, x, A_mul_x);
97  vector_t resid = vector_get_residual(i+1, A_mul_x, b);
98  printf("residual = %e\n", resid);
99 
100  vector_clear(i+1, x);
101  algo = Givens;
102  puts("Givens: ");
103 // matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
104 // i + 1, sub_A);
105  // Start Timer !!!
106  error_no = solve(i + 1, i + 1, sub_A, b, x, algo);
107  // Stop Timer !!!
108 
109  if (error_no > 0) {
110  printf("x = ");
111  vector_flex_print(i + 1, x, 3, 4);
112  puts("");
113  }
114 
115  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
116  i + 1, sub_A);
117  matrix_mul_vec(i+1, i+1, sub_A, x, A_mul_x);
118  resid = vector_get_residual(i+1, A_mul_x, b);
119  printf("residual = %e\n", resid);
120 
121  vector_clear(i+1, x);
122  algo = Householder;
123  puts("Householder: ");
124 // matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
125 // i + 1, sub_A);
126  // Start Timer !!!
127  error_no = solve(i + 1, i + 1, sub_A, b, x, algo);
128  // Stop Timer !!!
129 
130  if (error_no > 0) {
131  printf("x = ");
132  vector_flex_print(i + 1, x, 3, 4);
133  puts("");
134  }
135 
136  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
137  i + 1, sub_A);
138  matrix_mul_vec(i+1, i+1, sub_A, x, A_mul_x);
139  resid = vector_get_residual(i+1, A_mul_x, b);
140  printf("residual = %e\n", resid);
141 
142  vector_clear(i+1, x);
143  algo = Gauss;
144  puts("Gauss");
145 // matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
146 // i + 1, sub_A);
147  // Start Timer !!!
148  error_no = solve(i + 1, i + 1, sub_A, b, x, algo);
149  // Stop Timer !!!
150 
151  if (error_no > 0) {
152  printf("x = ");
153  vector_flex_print(i + 1, x, 3, 4);
154  puts("");
155  }
156 
157  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, i, i + 1,
158  i + 1, sub_A);
159  matrix_mul_vec(i+1, i+1, sub_A, x, A_mul_x);
160  resid = vector_get_residual(i+1, A_mul_x, b);
161  printf("residual = %e\n", resid);
162  } //for
163 
164  puts("Rectangular sub-matrices:");
165  for (uint8_t i = COL_MAX; i < ROW_MAX; i++) {
166  matrix_t sub_A[i + 1][COL_MAX];
167  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, COL_MAX - 1,
168  i + 1, COL_MAX, sub_A);
169  printf("(%u,%u) matrix: ", i + 1, COL_MAX);
170  matrix_flex_print(i + 1, COL_MAX, sub_A, 3, 4);
171  puts("");
172 
173  enum ALGORITHM algo = Moore_Penrose;
174  puts("Moore Penrose: ");
175  matrix_t x[COL_MAX];
176 
177  vector_clear(COL_MAX, x);
178  // Start Timer !!!
179  int8_t error_no = solve(i + 1, COL_MAX, sub_A, b, x, algo);
180  // Stop Timer !!!
181  if (error_no > 0) {
182  printf("x = ");
183  vector_flex_print(COL_MAX, x, 3, 4);
184  puts("");
185  }
186 
187  vector_clear(COL_MAX, x);
188  algo = Givens;
189  puts("Givens: ");
190  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, COL_MAX - 1,
191  i + 1, COL_MAX, sub_A);
192  // Start Timer !!!
193  error_no = solve(i + 1, COL_MAX, sub_A, b, x, algo);
194  // Stop Timer !!!
195  if (error_no > 0) {
196  printf("x = ");
197  vector_flex_print(COL_MAX, x, 3, 4);
198  puts("");
199  }
200 
201  vector_clear(COL_MAX, x);
202  algo = Householder;
203  puts("Householder: ");
204  matrix_part_copy(ROW_MAX, COL_MAX, copy_A, 0, i, 0, COL_MAX - 1,
205  i + 1, COL_MAX, sub_A);
206  // Start Timer !!!
207  error_no = solve(i + 1, COL_MAX, sub_A, b, x, algo);
208  // Stop Timer !!!
209  if (error_no > 0) {
210  printf("x = ");
211  vector_flex_print(COL_MAX, x, 3, 4);
212  puts("");
213  }
214  }//for
215 }
216 
217 void ieee_mobile_comp_test_qr_givens(void)
218 {
219  for (uint8_t i = 4; i < ROW_MAX; i++) {
220  matrix_t Q[i + 1][i + 1];
221  matrix_t red_Q[i + 1][COL_MAX];
222  matrix_t red_R[COL_MAX][COL_MAX];
223  matrix_part_copy(ROW_MAX, COL_MAX, A, 0, i, 0, COL_MAX, ROW_MAX,
224  COL_MAX, copy_A);
225  // printf("copy_A = ");
226  // matrix_flex_part_print(ROW_MAX, COL_MAX, copy_A, 0, i, 0, COL_MAX-1, 1, 4);
227  // puts("");
228  // Start Timer !!!
229  qr_givens_decomp(i + 1, COL_MAX, copy_A, COL_MAX, Q, false);
230  // Stop Timer !!!
231  printf("Decomposition of A_%d_5\n", i + 1);
232  qr_common_get_reduced_QR(i + 1, COL_MAX, Q, copy_A, red_Q,
233  red_R);
234  printf("Q = ");
235  matrix_flex_print(i + 1, COL_MAX, red_Q, 1, 4);
236  puts("");
237  printf("R = ");
238  matrix_flex_part_print(i, COL_MAX, red_R, 0, COL_MAX - 1, 0,
239  COL_MAX - 1, 1, 4);
240  puts("");
241  matrix_copy(ROW_MAX, COL_MAX, A, copy_A);
242  }
243 }
244 
245 
246 
247 
248 void ieee_mobile_comp_test_qr_householder(void)
249 {
250  for (uint8_t i = 4; i < ROW_MAX; i++) {
251  matrix_t Q[i + 1][i + 1];
252  matrix_t red_Q[i + 1][COL_MAX];
253  matrix_t red_R[COL_MAX][COL_MAX];
254  matrix_part_copy(ROW_MAX, COL_MAX, A, 0, i, 0, COL_MAX, ROW_MAX,
255  COL_MAX, copy_A);
256  // printf("copy_A = ");
257  // matrix_flex_part_print(ROW_MAX, COL_MAX, copy_A, 0, i, 0, COL_MAX-1, 1, 4);
258  // puts("");
259  // Start Timer !!!
260  qr_householder_decomp(i + 1, COL_MAX, copy_A, COL_MAX, Q, false);
261  // Stop Timer !!!
262  printf("Decomposition of A_%d_5\n", i + 1);
263  qr_common_get_reduced_QR(i + 1, COL_MAX, Q, copy_A, red_Q,
264  red_R);
265  printf("Q = ");
266  matrix_flex_print(i + 1, COL_MAX, red_Q, 1, 4);
267  puts("");
268  printf("R = ");
269  matrix_flex_part_print(i, COL_MAX, red_R, 0, COL_MAX - 1, 0,
270  COL_MAX - 1, 1, 4);
271  puts("");
272  matrix_copy(ROW_MAX, COL_MAX, A, copy_A);
273  }
274 }
275 
276 void ieee_mobile_comp_test_SVD(void)
277 {
278  //uint8_t m = ROW_MAX;
279  uint8_t n = COL_MAX;
280  matrix_dim_t dim;
281 
282  // svd_get_U_dim(m, n, &dim);
283  // matrix_t U2[dim.row_num][dim.col_num];
284  // matrix_t S2[n][n];
285  // matrix_t V2[n][n];
286  // uint8_t s2_length = svd_get_single_values_num(m, n);
287  // matrix_t s2[s2_length];
288  // svd(m, n, A, dim.row_num, dim.col_num, U2, S2, V2, s2_length, s2);
289  // printf("U2 = ");
290  // matrix_print(dim.row_num, dim.col_num, U2);
291  // puts("");
292  // printf("S2 = ");
293  // matrix_print(n, n, S2);
294  // puts("");
295  // printf("V2 = ");
296  // matrix_print(n, n, V2);
297  // puts("");
298 
299  for (uint8_t i = 4; i < ROW_MAX; i++) {
300  svd_get_U_dim(i + 1, COL_MAX, &dim);
301  matrix_t U[dim.row_num][dim.col_num];
302  matrix_t S[COL_MAX][COL_MAX];
303  matrix_t V[COL_MAX][COL_MAX];
304  uint8_t s_length = svd_get_single_values_num(i, COL_MAX);
305  matrix_t s[s_length];
306  matrix_part_copy(ROW_MAX, COL_MAX, A, 0, i, 0, COL_MAX, ROW_MAX,
307  COL_MAX, copy_A);
308  // printf("copy_A = ");
309  // matrix_flex_part_print(ROW_MAX, COL_MAX, copy_A, 0, i, 0,
310  // COL_MAX - 1, 1, 4);
311  // puts("");
312  svd(i + 1, COL_MAX, copy_A, dim.row_num, dim.col_num, U, S, V,
313  s_length,
314  s);
315  printf("U = ");
316  matrix_print(dim.row_num, dim.col_num, U);
317  puts("");
318  printf("S2 = ");
319  matrix_print(n, n, S);
320  puts("");
321  printf("V2 = ");
322  matrix_print(n, n, V);
323  puts("");
324 
325  }
326 }
327 
matrix_mul_vec
void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t vec[n], matrix_t dst_arr[m])
Compute the multiplication of a matrix with a column vector.
Definition: matrix.c:434
Givens
Givens algorithm.
Definition: pseudo_inverse.h:33
qr_householder_decomp
int8_t qr_householder_decomp(uint8_t m, uint8_t n, matrix_t A[][n], uint8_t q_col_num, matrix_t Q[][q_col_num], bool reduced)
Computes the QR decomposition of the matrix A by using the Householder algorithm.
Definition: qr_householder.c:33
vector_flex_print
void vector_flex_print(uint32_t length, vector_t arr[], uint8_t before_dot, uint8_t after_dot)
Display the values of the vector's elements.
Definition: vector.c:284
vector_get_residual
vector_t vector_get_residual(uint8_t length, vector_t a_vec[], vector_t b_vec[])
Compute the residual of two vectors.
Definition: vector.c:227
qr_householder.h
Householder algorithm for the QR-decomposition.
matrix_dim_t::row_num
uint8_t row_num
the row number
Definition: matrix.h:61
Moore_Penrose
Moore–Penrose algorithm.
Definition: pseudo_inverse.h:31
qr_common_get_reduced_QR
void qr_common_get_reduced_QR(uint8_t m, uint8_t n, matrix_t Q[m][m], matrix_t R[m][n], matrix_t red_Q[m][n], matrix_t red_R[n][n])
Compute the reduced form of the QR-decomposition algorithm.
Definition: qr_common.c:59
svd
void svd(uint8_t m, uint8_t n, matrix_t A[m][n], uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n], matrix_t S[u_n][n], matrix_t V[n][n], uint8_t sing_vec_length, matrix_t singl_values_vec[])
Compute the Singular-Value Decomposition (SVD) of a matrix.
Definition: svd.c:119
matrix_part_mul
void matrix_part_mul(uint8_t a_col_num_max, matrix_t a_matrix[][a_col_num_max], uint8_t b_col_num_max, matrix_t b_matrix[][b_col_num_max], uint8_t a_start_row_ind, uint8_t a_end_row_ind, uint8_t a_start_col_ind, uint8_t a_end_col_ind, uint8_t b_start_row_ind, uint8_t b_end_row_ind, uint8_t b_start_col_ind, uint8_t b_end_col_ind, uint8_t dest_col_size, matrix_t dest_matrix[][dest_col_size])
Compute the partial multiplication of two matrices.
Definition: matrix.c:391
vector_t
#define vector_t
Define the data type of the vector elements.
Definition: vector.h:33
solve.h
Enables to solve systems of linear equations Ax = b for x.
matrix.h
Matrix computations.
matrix_flex_print
void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t before_dec, uint8_t after_dec)
Display the values of the matrix elements.
Definition: matrix.c:220
Householder
Householder algorithm.
Definition: pseudo_inverse.h:32
matrix_dim_t
A structure to define the row and column number of a matrix.
Definition: matrix.h:60
matrix_flex_part_print
void matrix_flex_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t start_row_ind, uint8_t end_row_ind, uint8_t start_col_ind, uint8_t end_col_ind, uint8_t before_dot, uint8_t after_dot)
Display the values of the sub-matrix elements.
Definition: matrix.c:247
vector_clear
void vector_clear(uint8_t size, vector_t arr[])
Clear all the elements of the vector.
Definition: vector.c:32
ALGORITHM
ALGORITHM
Possible algorithms to compute the pseudo-inverse matrix.
Definition: pseudo_inverse.h:30
solve
int8_t solve(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m], matrix_t x_sol[n], enum ALGORITHM algo)
Solve an (m n) linear system Ax = b by using the Moore–Penrose, Householder, or the Givens algorithm...
Definition: solve.c:35
qr_common.h
Common definitions and implementations for the QR-decomposition. Provide necessary methods to constru...
svd.h
Algorithm for the Singular Value Decomposition (SVD).
Gauss
Gaussian elimination with pivoting algorithm.
Definition: pseudo_inverse.h:34
svd_get_single_values_num
uint8_t svd_get_single_values_num(uint8_t m, uint8_t n)
Calculate the number of the singular values.
Definition: svd.c:114
matrix_copy
void matrix_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], matrix_t dest_matrix[m][n])
Copy the elements of a matrix to another matrix.
Definition: matrix.c:83
qr_givens_decomp
int8_t qr_givens_decomp(uint8_t m, uint8_t n, matrix_t A[][n], uint8_t q_col_num, matrix_t Q[][q_col_num], bool reduced)
Computes the QR decomposition of the matrix A by using the Givens algorithm.
Definition: qr_givens.c:33
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
vector.h
Vector computations.
matrix_print
void matrix_print(uint8_t m, uint8_t n, matrix_t matrix[m][n])
Display the values of the matrix elements.
Definition: matrix.c:141
matrix_part_copy
void matrix_part_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], uint8_t start_row_ind, uint8_t end_row_ind, uint8_t start_col_ind, uint8_t end_col_ind, uint8_t dest_row_num, uint8_t dest_col_num, matrix_t dest_matrix[][dest_col_num])
Copy a part of a matrix to another matrix or sub-matrix.
Definition: matrix.c:89
qr_givens.h
Givens algorithm for the QR-decomposition. Provide necessary methods to construct Q- and R- matrices ...
svd_get_U_dim
void svd_get_U_dim(uint8_t m, uint8_t n, matrix_dim_t *u_dim)
Calculate the dimension of the matrix U.
Definition: svd.c:96
matrix_dim_t::col_num
uint8_t col_num
the column number
Definition: matrix.h:62