RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
matrix_test.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
3  * 2020 Freie Universität Berlin
4  *
5  * This file is subject to the terms and conditions of the GNU Lesser General
6  * Public License v2.1. See the file LICENSE in the top level directory for more
7  * details.
8  */
9 
23 #include <stdio.h>
24 #include <stdint.h>
25 #include <inttypes.h>
26 #include <math.h>
27 
28 #include "utils.h"
29 #include "matrix.h"
30 
31 void matrix_test(void)
32 {
33  puts("############ Basic Matrix Algebra ###############");
34  uint8_t m = 3;
35  uint8_t n = 4;
36  matrix_t diag_elem = M_PI;
37 
38  matrix_t matrix1[3][4] = { { 0.2785, 0.9649, 0.9572, 0.1419 },
39  { 0.5469, 0.1576, 0.4854, 0.4218 },
40  { 0.9575, 0.9706, 0.8003, 0.9157 },
41  };
42 
43  matrix_t matrix2[3][4] = { { 0.7922, 0.0357, 0.6787, 0.3922 },
44  { 0.9595, 0.8491, 0.7577, 0.6555 },
45  { 0.6557, 0.9340, 0.7431, 0.1712 },
46  };
47  matrix_t res_matrix[m][n];
48  matrix_t trans_matrix[n][m];
49  matrix_t res_mul_matrix[m][m];
50  matrix_t diag_vec[3] = { 1.1, 2.3, 3.5 };
51 
52  matrix_add(m, n, matrix1, matrix2, res_matrix);
53  printf("matrix1 + matrix2 = ");
54  matrix_flex_print(m, n, res_matrix, 7, 4);
55 
56  matrix_sub(m, n, matrix1, matrix2, res_matrix);
57  printf("matrix1 - matrix2 = ");
58  matrix_flex_print(m, n, res_matrix, 7, 4);
59 
60  matrix_transpose(m, n, matrix2, trans_matrix);
61  printf("trans(matrix2) = ");
62  matrix_flex_print(n, m, trans_matrix, 7, 4);
63 
64  matrix_mul(m, n, matrix1, n, m, trans_matrix, res_mul_matrix);
65  printf("matrix1 x matrix2 = ");
66  matrix_flex_print(m, m, res_mul_matrix, 7, 4);
67 
68  matrix_get_diag_mat(m, n, diag_elem, res_matrix);
69  printf("dig_matrix = ");
70  matrix_flex_print(m, n, res_matrix, 7, 4);
71 
72  matrix_get_diag_mat_new(m, n, res_matrix, m, diag_vec);
73  printf("dig_matrix = ");
74  matrix_flex_print(m, n, res_matrix, 7, 4);
75 
76  matrix_t matrix3[5][5] = {
77  { 0.8147, 0.0975, 0.1576, 0.1419, 0.6557 },
78  { 0.9058, 0.2785, 0.9706, 0.4218, 0.0357 },
79  { 0.1270, 0.5469, 0.9572, 0.9157, 0.8491 },
80  { 0.9595, 0.8491, 0.7577, 0.6555, 0.6787 },
81  { 0.6557, 0.9340, 0.7431, 0.1712, 0.9595 },
82  };
83 
84  matrix_swap_rows(5, matrix3, 3, 1);
85  printf("swapped matrix3 = ");
86  matrix_flex_print(5, 5, matrix3, 7, 4);
87 
88  puts(" *** Test Partial Matrix Multiplication *** ");
89  matrix_t A[11][7] =
90  {
91  { 0.8147, 0.9706, 0.8491, 0.0462, 0.1869, 0.4984, 0.5472 },
92  { 0.9058, 0.9572, 0.9340, 0.0971, 0.4898, 0.9597, 0.1386 },
93  { 0.1270, 0.4854, 0.6787, 0.8235, 0.4456, 0.3404, 0.1493 },
94  { 0.9134, 0.8003, 0.7577, 0.6948, 0.6463, 0.5853, 0.2575 },
95  { 0.6324, 0.1419, 0.7431, 0.3171, 0.7094, 0.2238, 0.8407 },
96  { 0.0975, 0.4218, 0.3922, 0.9502, 0.7547, 0.7513, 0.2543 },
97  { 0.2785, 0.9157, 0.6555, 0.0344, 0.2760, 0.2551, 0.8143 },
98  { 0.5469, 0.7922, 0.1712, 0.4387, 0.6797, 0.5060, 0.2435 },
99  { 0.9575, 0.9595, 0.7060, 0.3816, 0.6551, 0.6991, 0.9293 },
100  { 0.9649, 0.6557, 0.0318, 0.7655, 0.1626, 0.8909, 0.3500 },
101  { 0.1576, 0.0357, 0.2769, 0.7952, 0.1190, 0.9593, 0.1966 }
102  };
103 
104  matrix_t B[7][5] = {
105  { 0.2511, 0.9172, 0.0540, 0.0119, 0.6020 },
106  { 0.6160, 0.2858, 0.5308, 0.3371, 0.2630 },
107  { 0.4733, 0.7572, 0.7792, 0.1622, 0.6541 },
108  { 0.3517, 0.7537, 0.9340, 0.7943, 0.6892 },
109  { 0.8308, 0.3804, 0.1299, 0.3112, 0.7482 },
110  { 0.5853, 0.5678, 0.5688, 0.5285, 0.4505 },
111  { 0.5497, 0.0759, 0.4694, 0.1656, 0.0838 }
112  };
113 
114  uint8_t l;
115  uint8_t a_row_begin, a_row_end, a_col_begin, a_col_end;
116  uint8_t b_row_begin, b_row_end, b_col_begin, b_col_end;
117  m = 11;
118  n = 7;
119  l = 5;
120  matrix_t C[5][3];
121 
122  a_row_begin = 3, a_row_end = 7;
123  a_col_begin = 2, a_col_end = 4;
124  b_row_begin = 3, b_row_end = 5;
125  b_col_begin = 1, b_col_end = 3;
126  matrix_part_mul(n, A, l, B,
127  a_row_begin, a_row_end, a_col_begin, a_col_end,
128  b_row_begin, b_row_end, b_col_begin, b_col_end,
129  3, C);
130 
131  printf("part_A = ");
132  matrix_flex_part_print(m, n, A, a_row_begin, a_row_end, a_col_begin,
133  a_col_end, 7, 4);
134  printf("part_B = ");
135  matrix_flex_part_print(n, l, B, b_row_begin, b_row_end, b_col_begin,
136  b_col_end, 7, 4);
137 
138  printf("C = ");
139  matrix_flex_print(5, 3, C, 7, 4);
140 
141  puts(" *** Test The 2-norm of a Matrix, *** ");
142 
143  double two_norm = matrix_get_two_norm(m, n, A);
144  printf("A_two_norm = %7.4f\n", two_norm);
145 
146  matrix_t D[7][5] = {
147  { -0.2743, -4.9913, 0.3375, -6.5684, -6.9048 },
148  { -4.9657, -2.0716, -4.7133, -2.7536, -4.3030 },
149  { -0.4857, -3.2137, -0.9424, -0.7667, -5.7025 },
150  { -5.0518, -4.1867, -0.9702, 0.4721, -0.6457 },
151  { 0.4341, -0.3534, -3.9564, -5.9608, -4.5103 },
152  { -4.2001, -2.3179, -2.4574, -2.4494, -2.7717 },
153  { -5.4272, -2.6022, -6.3932, -3.2449, -5.6748 }
154  };
155  m = 7;
156  n = 5;
157  two_norm = matrix_get_two_norm(m, n, D);
158  printf("D_two_norm = %7.4f\n", two_norm);
159 
160  double frob_norm = matrix_get_frob_norm(m, n, D);
161  printf("D_frob_norm = %7.4f\n", frob_norm);
162 }
163 
165 {
166 
167  puts(
168  "*********** Test the Inverse of Upper Triangular Matrices ***********");
169  matrix_t U[11][7] =
170  {
171  { 0.8147, 0.9706, 0.8491, 0.0462, 0.1869, 0.4984, 0.5472 },
172  { 0.0000, 0.9572, 0.9340, 0.0971, 0.4898, 0.9597, 0.1386 },
173  { 0.0000, 0.0000, 0.6787, 0.8235, 0.4456, 0.3404, 0.1493 },
174  { 0.0000, 0.0000, 0.0000, 0.6948, 0.6463, 0.5853, 0.2575 },
175  { 0.0000, 0.0000, 0.0000, 0.0000, 0.7094, 0.2238, 0.8407 },
176  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.7513, 0.2543 },
177  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.8143 },
178  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 },
179  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 },
180  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 },
181  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 }
182  };
183  uint8_t m, n;
184 
185  puts("Square Matrix:");
186  m = 7;
187  n = 7;
188  printf("U1 = ");
189  matrix_flex_print(m, n, U, 7, 4);
190  matrix_t inv_U1[n][m];
191  matrix_get_inv_upp_triang(m, n, U, inv_U1);
192  printf("inv_U1 = ");
193  matrix_flex_print(n, m, inv_U1, 7, 4);
194 
195  puts("Rectangular Matrix:");
196  m = 11;
197  n = 7;
198  printf("U2 = ");
199  matrix_flex_print(m, n, U, 7, 4);
200  matrix_t inv_U2[n][m];
201  matrix_get_inv_upp_triang(m, n, U, inv_U2);
202  printf("inv_U2 = ");
203  matrix_flex_print(n, m, inv_U2, 7, 4);
204 
205  puts(
206  "*********** Test the Inverse of Lower Triangular Matrices ***********");
207  matrix_t L[11][7] =
208  {
209  { 0.8147, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 },
210  { 0.9058, 0.9572, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 },
211  { 0.1270, 0.4854, 0.6787, 0.0000, 0.0000, 0.0000, 0.0000 },
212  { 0.9134, 0.8003, 0.7577, 0.6948, 0.0000, 0.0000, 0.0000 },
213  { 0.6324, 0.1419, 0.7431, 0.3171, 0.7094, 0.0000, 0.0000 },
214  { 0.0975, 0.4218, 0.3922, 0.9502, 0.7547, 0.7513, 0.0000 },
215  { 0.2785, 0.9157, 0.6555, 0.0344, 0.2760, 0.2551, 0.8143 },
216  { 0.5469, 0.7922, 0.1712, 0.4387, 0.6797, 0.5060, 0.2435 },
217  { 0.9575, 0.9595, 0.7060, 0.3816, 0.6551, 0.6991, 0.9293 },
218  { 0.9649, 0.6557, 0.0318, 0.7655, 0.1626, 0.8909, 0.3500 },
219  { 0.1576, 0.0357, 0.2769, 0.7952, 0.1190, 0.9593, 0.1966 }
220  };
221  m = 7;
222  n = 7;
223  puts("Square Matrix:");
224  printf("L1 = ");
225  matrix_flex_print(m, n, L, 7, 4);
226  matrix_t inv_L1[n][m];
227  matrix_get_inv_low_triang(m, n, L, inv_L1);
228  printf("inv_L1 = ");
229  matrix_flex_print(n, m, inv_L1, 7, 4);
230 
231  puts("Rectangular Matrix:");
232  m = 11;
233  n = 7;
234  printf("L2 = ");
235  matrix_flex_print(m, n, L, 7, 4);
236  matrix_t inv_L2[n][m];
237  matrix_get_inv_low_triang(m, n, L, inv_L2);
238  printf("inv_L2 = ");
239  matrix_flex_print(n, m, inv_L2, 7, 4);
240 }
241 
243 {
244  matrix_t A[10][5] = {
245  { 0.8147, 0.1576, 0.6557, 0.7060, 0.4387 },
246  { 0.9058, 0.9706, 0.0357, 0.0318, 0.3816 },
247  { 0.1270, 0.9572, 0.8491, 0.2769, 0.7655 },
248  { 0.9134, 0.4854, 0.9340, 0.0462, 0.7952 },
249  { 0.6324, 0.8003, 0.6787, 0.0971, 0.1869 },
250  { 0.0975, 0.1419, 0.7577, 0.8235, 0.4898 },
251  { 0.2785, 0.4218, 0.7431, 0.6948, 0.4456 },
252  { 0.5469, 0.9157, 0.3922, 0.3171, 0.6463 },
253  { 0.9575, 0.7922, 0.6555, 0.9502, 0.7094 },
254  { 0.9649, 0.9595, 0.1712, 0.0344, 0.7547 }
255  };
256 
257  uint8_t m = 10, n = 5;
258  matrix_t copy_A[m][n];
259 
260  matrix_copy(m, n, A, copy_A);
261  matrix_t B[m][n];
262  matrix_t C[n][m];
263  matrix_t D[n][m];
264 
265  puts("++++++++++ upper triangular ++++++++++");
266  matrix_transpose(m, n, A, C);
267  matrix_get_upp_triang(m, n, A, B);
268  matrix_get_upp_triang(m, n, A, A);
269  matrix_get_upp_triang(n, m, C, D);
270 
271  printf("A = ");
272  matrix_flex_print(m, n, A, 7, 4);
273  printf("B = ");
274  matrix_flex_print(m, n, B, 7, 4);
275  printf("C = ");
276  matrix_flex_print(n, m, C, 7, 4);
277  printf("D = ");
278  matrix_flex_print(n, m, D, 7, 4);
279 
280  puts("++++++++++ lower triangular ++++++++++");
281  matrix_copy(m, n, copy_A, A);
282  printf("A = ");
283  matrix_flex_print(m, n, A, 7, 4);
284  matrix_transpose(m, n, A, C);
285  matrix_get_low_triang(m, n, A, B);
286  matrix_get_low_triang(m, n, A, A);
287  matrix_get_low_triang(n, m, C, D);
288 
289  printf("A = ");
290  matrix_flex_print(m, n, A, 7, 4);
291  printf("B = ");
292  matrix_flex_print(m, n, B, 7, 4);
293  printf("C = ");
294  matrix_flex_print(n, m, C, 7, 4);
295  printf("D = ");
296  matrix_flex_print(n, m, D, 7, 4);
297 
298 }
matrix_swap_rows
void matrix_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j)
Swaps two rows of a matrix.
Definition: matrix.c:725
matrix_sub
void matrix_sub(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n], matrix_t A_minus_B[m][n])
Compute the subtraction of two matrices.
Definition: matrix.c:329
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
matrix_transpose
void matrix_transpose(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], matrix_t dest_matrix[n][m])
Computes the transpose of a matrix.
Definition: matrix.c:56
matrix_add
void matrix_add(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n], matrix_t A_plus_B[m][n])
Compute the addition of two matrices.
Definition: matrix.c:341
triangular_matrices_test
void triangular_matrices_test(void)
Examples of triangular matrices.
Definition: matrix_test.c:242
matrix_test
void matrix_test(void)
Test some matrix operations.
Definition: matrix_test.c:31
matrix_get_upp_triang
void matrix_get_upp_triang(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t tr_up_A[][n])
Gets the upper triangular part of a matrix.
Definition: matrix.c:841
inv_triangular_matrices_test
void inv_triangular_matrices_test(void)
Examples of the inverse of triangular matrices.
Definition: matrix_test.c:164
matrix.h
Matrix computations.
matrix_get_diag_mat
void matrix_get_diag_mat(uint8_t m, uint8_t n, matrix_t value, matrix_t diag_matrix[m][n])
Create a diagonal matrix with a specified value.
Definition: matrix.c:594
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
utils.h
Utilities for linear algebra.
matrix_get_low_triang
void matrix_get_low_triang(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t tr_low_A[][n])
Gets the lower triangular part of a matrix.
Definition: matrix.c:868
M_PI
#define M_PI
Definition: matrix.h:54
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
matrix_get_inv_low_triang
void matrix_get_inv_low_triang(uint8_t m, uint8_t n, matrix_t L[][n], matrix_t inv_L[][m])
Computes the inverse a lower triangular matrix.
Definition: matrix.c:817
matrix_get_two_norm
double matrix_get_two_norm(uint8_t m, uint8_t n, matrix_t A[][n])
Get the 2-norm of a matrix that is equal to the largest singular value.
Definition: matrix.c:752
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
matrix_get_inv_upp_triang
void matrix_get_inv_upp_triang(uint8_t m, uint8_t n, matrix_t U[][n], matrix_t inv_U[][m])
Computes the inverse an upper triangular matrix.
Definition: matrix.c:795
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
matrix_get_frob_norm
double matrix_get_frob_norm(uint8_t m, uint8_t n, matrix_t A[][n])
Get the Frobenius norm of a matrix.
Definition: matrix.c:782
matrix_mul
void matrix_mul(uint8_t a_line_num, uint8_t a_col_num, matrix_t a_matrix[a_line_num][a_col_num], uint8_t b_line_num, uint8_t b_col_num, matrix_t b_matrix[b_line_num][b_col_num], matrix_t dest_matrix[a_line_num][b_col_num])
Compute the multiplication of two matrices.
Definition: matrix.c:363
matrix_get_diag_mat_new
void matrix_get_diag_mat_new(uint8_t m, uint8_t n, matrix_t diag_matrix[m][n], uint8_t length, matrix_t vec[])
Set all the diagonal elements of a matrix with values that are saved in a vector.
Definition: matrix.c:581